<!DOCTYPE HTML PUBLIC "-//SoftQuad Software//DTD HoTMetaL PRO 5.0::19981217::extensions to HTML 4.0//EN" "hmpro5.dtd">
 
<HTML>
 
<HEAD>
<meta http-equiv="Content-Language" content="en-nz">
<meta name="GENERATOR" content="Microsoft FrontPage 6.0">
<meta name="ProgId" content="FrontPage.Editor.Document">

<TITLE>Newmat11 documentation</TITLE>
<link rel="stylesheet" type="text/css" href="rbd.css">
<link REL="SHORTCUT ICON" HREF="favicon.ico">
</HEAD>

<BODY>
<H1><A NAME="top"></A>Documentation for newmat11, a matrix library in C++</H1>
<P CLASS="small"><A HREF="#intro">next</A> - <A HREF="#top">skip</A> -
<A HREF="ol_doc.htm">up</A> - <A HREF="#top">start</A><BR>
<A HREF="ol_doc.htm">return to online documentation page</A></P>
<P><B>Copyright (C)  2008: R B Davies</B></P>
<P><i>20 November, 2008</i>.</P>
<TABLE WIDTH=100%>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#intro">1.
Introduction</A><BR>
<A HREF="#starting">2. Getting started</A><BR>
<A HREF="#refer">3. Reference manual</A><br>
<A HREF="#error">4.
Error messages</A></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> 
<A HREF="#design">5. Design of the library</A><br>
<a href="#function">6. Function summary</a><br>
<a href="#changes">7. Change History</a><br>
<a href="#problem">8. Problem report form</a></TD>
</TR>
</TABLE>
<p>&nbsp;</p>
<P>This is the <I>how to use</I> documentation for <I>newmat</I> plus some
background information on its design. </P>
<P>There is additional support material on my <A HREF="#where">web site</A>. 
</P>
<P CLASS="small">Navigation:&nbsp; This page is arranged in sections,
sub-sections and sub-sub-sections; four cross-references are given at the top
of these. <I>Next</I> takes you through the sections, sub-sections and
sub-sub-sections in order. <I>Skip</I> goes to the next section, sub-section or
sub-sub-section at the same level in the hierarchy as the section, sub-section
or sub-sub-section that you are currently reading. <I>Up</I> takes you up one
level in the hierarchy and <I>start</I> gets you back here.</P>
<P><B>Please read the sections on <A HREF="#custom">customising</A> and
<a href="#compiler">compilers</a> before
attempting to compile <i>newmat</i>.</B> </P>
<H2><A NAME="intro"></A>1. Introduction</H2>
<P CLASS="small"><A HREF="#use">next</A> - <A HREF="#starting">skip</A> -
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH=100%>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#use">1.1
Conditions of use</A><BR>
<A HREF="#descript">1.2 Description</A><BR>
<A HREF="#which">1.3 Is this the library for you?</A><BR>
<A HREF="#other">1.4 Other matrix libraries</A></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#where">1.5
Where to find this library</A><BR>
<A HREF="#author">1.6 How to contact the author</A><BR>
<A HREF="#sources">1.7 References</A></TD>
</TR>
</TABLE>
<H2><A NAME="use"></A>1.1 Conditions of use</H2>
<P CLASS="small"><A HREF="#descript">next</A> - <A HREF="#descript">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<hr>
<p>I place no restrictions on the use of <i>newmat</i> except that I take no liability for any problems that may arise from its use,
distribution or other dealings with it.<p>
<p>
You can use it in your commercial projects (as well as your non-commercial 
projects).<p>
<p>
You can make and distribute modified or merged versions. You can
include parts of it in your own software.<p>
<p>
If you distribute modified or merged versions, please make it clear
which parts are mine and which parts are modified.<p>
<p>
For a substantially modified version, simply note that it is, in
part, derived from my software. A comment in the code will be
sufficient.<p>
<p>
The software is provided <i>as is</i>, without warranty of any kind.<p>
<p>
Please understand that there may still be bugs and errors. Use at
your own risk. I (Robert Davies) take no responsibility for any errors
or omissions in this package or for any misfortune that may befall you
or others as a result of your use, distribution or other dealings with it.<HR>
<P>Please report bugs to me at <B>robert at statsresearch.co.nz</B>&nbsp;&nbsp; 
[replace <b>at</b> by you-know-what-character in the email address].</P>
<P>When reporting a bug please tell me which C++ compiler you are using, and
what version. Also give me details of your computer. And tell me which version
of <I>newmat</I> (e.g. newmat03 or newmat04) you are using and its date. Note any changes
you have made to my code. If at all possible give me a piece of code
illustrating the bug. See the <A HREF="#problem">problem report form</A>. </P>
<P><I>Please do report bugs to me.</I> </P>
<H2><A NAME="descript"></A>1.2 General description
</H2>
<P CLASS="small"><A HREF="#which">next</A> - <A HREF="#which">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<P>The package is intended for scientists and engineers who need to manipulate
a variety of types of matrices using standard matrix operations. Emphasis is on
the kind of operations needed in statistical calculations such as least
squares, linear equation solve and eigenvalues. </P>
<P>It supports matrix types </P>
<BLOCKQUOTE>
  <TABLE WIDTH="80%" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111" cellpadding="0">
<TR>
<TD>Matrix</TD>
<TD>rectangular matrix</TD>
</TR>
<TR>
<TD>SquareMatrix</TD>
<TD>square matrix</TD>
</TR>
<TR>
<TD>nricMatrix</TD>
<TD>for use with <I>Numerical Recipes in C</I> programs</TD>
</TR>
<TR>
<TD>UpperTriangularMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>LowerTriangularMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>DiagonalMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>SymmetricMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>BandMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>UpperBandMatrix</TD>
<TD>upper triangular band matrix</TD>
</TR>
<TR>
<TD>LowerBandMatrix</TD>
<TD>lower triangular band matrix</TD>
</TR>
<TR>
<TD>SymmetricBandMatrix</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>RowVector</TD>
<TD>derived from Matrix</TD>
</TR>
<TR>
<TD>ColumnVector</TD>
<TD>derived from Matrix</TD>
</TR>
<TR>
<TD>IdentityMatrix</TD>
<TD>diagonal matrix, elements have same value</TD>
</TR>
</TABLE>
</BLOCKQUOTE>
<P>Only one element type (float or double) is supported. </P>
<P>The package includes the operations <TT>*</TT>, <TT>+</TT>, <TT>-</TT>,
Kronecker product, Schur product, concatenation, inverse, transpose, conversion between types, submatrix,
determinant, Cholesky decomposition, QR decomposition, singular value
decomposition, eigenvalues of a symmetric matrix, sorting, fast Fourier
transform, printing and an interface with <I>Numerical Recipes in C</I>. </P>
<P>It is intended for matrices in the range 10 x 10 to the maximum size your
machine will accommodate in a single array. The number of elements in an array
cannot exceed the maximum size of an <I>int</I>. The package will work for very
small matrices but becomes rather inefficient. Some of the factorisation
functions are not (yet) optimised for paged memory and so become inefficient
when used with very large matrices. </P>
<P>A <I>lazy evaluation</I> approach to evaluating matrix expressions is used
to improve efficiency and reduce the use of temporary storage. </P>
<P>I have tested versions of the package on variety of compilers and platforms
including  Borland, Gnu, Microsoft and Sun. For more details
see the section on <A HREF="#compiler">compilers</A>. </P>
<H2><A NAME="which"></A>1.3 Is this the library for
you?</H2>
<P CLASS="small"><A HREF="#other">next</A> - <A HREF="#other">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<P>Do you</P>
<UL>
<LI>understand <TT>*</TT> to mean matrix multiply and not element by element
multiply </LI>
<LI>need matrix operators such as <TT>*</TT> and <TT>+</TT> defined as
operators so you can write things like <TT> X = A * (B + C);</TT></LI>
<LI>need a variety of types of matrices (but not sparse)</LI>
<LI>need only one element type (float or double)</LI>
<LI>work with matrices in the range 10 x 10 up to what can be stored in memory 
</LI>
<LI>tolerate a moderately large but not huge package</LI>
<LI>need high quality but not necessarily the latest numerical methods. </LI>
</UL>
<P>Then <I>newmat</I> may be the right matrix library for you. </P>
<H2><A NAME="other"></A>1.4 Other matrix libraries
</H2>
<P CLASS="small"><A HREF="#where">next</A> - <A HREF="#where">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<P>For details of other C++ matrix libraries look at
<a href="http://www.robertnz.net/cpp_site.html">http://www.robertnz.net/cpp_site.html</a>.
Look at the section <I>lists of libraries</I> which gives the locations of
several very comprehensive lists of matrix and other C++ libraries and at the
section <I>source code</I>. Or just search on <a href="http://www.google.com">Google</a>.</P>
<H2><A NAME="where"></A>1.5 Where to find this
library</H2>
<P CLASS="small"><A HREF="#author">next</A> - <A HREF="#author">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A> </P>
<UL>
<LI><a href="http://www.robertnz.net">http://www.robertnz.net</a></LI>
</UL>
<H2><A NAME="author"></A>1.6 How to contact the
author</H2>
<P CLASS="small"><a href="#sources">next</a> - <a href="#sources">skip</a> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<PRE>   Robert Davies
   16 Gloucester Street
   Wilton
   Wellington
   New Zealand

   <I>Internet:</I> robert <b>at</b> statsresearch.co.nz</PRE>

[replace <b>at</b> by you-know-what.]


<H2><A NAME="sources"></A>1.7 References</H2>
<P CLASS="small"><A HREF="#starting">next</A> - <A HREF="#starting">skip</A> -
<A HREF="#intro">up</A> - <A HREF="#top">start</A></P>
<UL>
<LI>The matrix LU decomposition is from Golub, G.H. &amp; Van Loan, C.F.
(1996), <I>Matrix Computations</I>, published by Johns Hopkins University
Press. </LI>
<LI>Part of the matrix inverse/solve routine is adapted from Press, Flannery,
Teukolsky, Vetterling (1988), <I>Numerical Recipes in C</I>, published by the
Cambridge University Press. </LI>
<LI>Many of the advanced matrix routines are adapted from routines in Wilkinson
and Reinsch (1971), <I>Handbook for Automatic Computation, Vol II, Linear
Algebra</I> published by Springer Verlag. </LI>
<LI>The fast Fourier transform is adapted from Carl de Boor (1980), <I>Siam J
Sci Stat Comput</I>, pp173-8 and the fast trigonometric transforms from Charles
Van Loan (1992) in <I>Computational frameworks for the fast Fourier
transform</I> published by SIAM. </LI>
<LI>The sort function is derived from Sedgewick, Robert (1992), <I>Algorithms
in C++</I> published by Addison Wesley. </LI>
</UL>
<P>For references about <I>Newmat</I> see </P>
<UL>
<LI>Davies, R.B. (1994) Writing a matrix package in C++. In OON-SKI'94: The
second annual object-oriented numerics conference, pp 207-213. Rogue Wave
Software, Corvallis. </LI>
<LI>Eddelbuttel, Dirk (1996) Object-oriented econometrics: matrix programming
in C++ using GCC and Newmat. Journal of Applied Econometrics, Vol 11, No 2, pp
199-209. </LI>
</UL>
<H2><A NAME="starting"></A>2. Getting started</H2>
<P CLASS="small"><A HREF="#overview">next</A> - <A HREF="#refer">skip</A> -
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH="100%">
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
HREF="#overview">2.1 Overview</A><BR>
<A HREF="#make">2.2 Make files</A><BR>
<A HREF="#custom">2.3 Customising</A><BR>
<A HREF="#compiler">2.4 Compilers</A><BR>
<a HREF="#update">2.5 Updating from previous versions</a><br>
<a HREF="#except_1">2.6 Catching exceptions</a></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
HREF="#example">2.7 Examples</A><BR>
<A HREF="#testing">2.8 Testing</A><BR>
<A HREF="#bugs">2.9 Bugs</A><BR>
<a href="#problemareas">2.10 Problem areas</a><BR>
<A HREF="#files">2.11 Files in newmat11</A><BR>
&nbsp;</TD>
</TR>
</TABLE>
<H2><A NAME="overview"></A>2.1 Overview</H2>
<P CLASS="small"><A HREF="#make">next</A> - <A HREF="#make">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>I use .h as the suffix of definition files and .cpp as the suffix of C++
source files. </P>
<P>You will need to compile all the *.cpp files listed as program files in the
<A HREF="#files">files section</A> to get the complete package. Ideally you
should store the resulting object files as a library. The tmt*.cpp files are
used for <A HREF="#testing">testing</A>, example.cpp is an <A
HREF="#example">example</A> and sl_ex.cpp, nl_ex.cpp and garch.cpp are examples
of the <A HREF="#nonlin">non-linear</A> solve and optimisation routines. A
demonstration and test of the exception mechanism is in test_exc.cpp. The files 
nm_ex1.cpp, nm_ex2.cpp and nm_ex3.cpp contain more simple examples. </P>
<P>I include a number of <I>make</I> files for compiling the example and the
test package. See the section on <A HREF="#make">make files</A> for details.
Alternatively, with the PC compilers, its pretty quick just to load all the files in the
interactive environments by pointing and clicking. </P>
<P>Use the large or win32 console model when you are using a PC. Do not
<I>outline</I> inline functions. You may need to increase the stack size on 
older operating systems or compilers. </P>
<P>Your source files that access the newmat will need to #include one or more
of the following files. </P>
<TABLE WIDTH="100%" CELLPADDING="3" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111">
<TR>
<TD VALIGN="TOP">include.h</TD>
<TD>to access just the compiler options</TD>
</TR>
<TR>
<TD VALIGN="TOP">newmat.h</TD>
<TD>to access just the main matrix library (includes include.h)</TD>
</TR>
<TR>
<TD VALIGN="TOP">newmatap.h</TD>
<TD>to access the advanced matrix routines such as Cholesky decomposition, QR
triangularisation etc (includes newmat.h)</TD>
</TR>
<TR>
<TD VALIGN="TOP">newmatio.h</TD>
<TD>to access the <A HREF="#output">output</A> routines (includes newmat.h) You
can use this only with compilers that support the standard input/output
routines including manipulators (all recent compilers)</TD>
</TR>
<TR>
<TD VALIGN="TOP">newmatnl.h</TD>
<TD>to access the non-linear optimisation routines (includes newmat.h)</TD>
</TR>
</TABLE>
<P>See the section on <A HREF="#custom">customising</A> to see how to edit
include.h for your environment and the section on <A
HREF="#compiler">compilers</A> for any special problems with the compiler you
are using.</P>
<H2><A NAME="make"></A>2.2 Make files</H2>
<P CLASS="small"><A HREF="#custom">next</A> - <A HREF="#custom">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>I have included <I>make</I> files for compiling my test and example programs for 
various versions of CC, Borland, Microsoft, Open Watcom, Intel and Gnu compilers. You 
can generate make files for a number of other compilers with my
<a href="genmake.htm">genmake</a> utility. <i>Make</i> files provide
a way of compiling your programs  without using the IDE that comes with
PC compilers. See the <A HREF="#files">files section</A> for details.</P>
<H3>PC</H3>
<P>I include  make files for various versions of Borland, Microsoft and Intel 
compilers. With the 
Borland compiler you will need to edit it to show where 
you have stored your Borland compiler. For make files for other compilers use my
<a href="genmake.htm">genmake</a> utility. To compile my test and example 
programs use Borland 5.5 (Builder 5) use</P>
<pre>   make -f nm_b55.mak</pre>
<p>or with Borland 5.6 (Builder 6) use</p>
<pre>   make -f nm_b56.mak</pre>
<P>or with Microsoft VC++ 6 or 7 use</P>
<pre>   nmake -f nm_m6.mak</pre>
<P>There are some more notes in the
<a href="genmake.htm">genmake</a> documentation about using these make files.</P>
<H3>Unix</H3>
<P>The <I>make</I> file for the Unix CC compilers link a .cxx file to each .cpp
file since some of these compilers do not recognise .cpp as a legitimate
extension for a C++ file. I suggest you delete this part of the <I>make</I>
file and, if necessary, rename the .cpp files to something your compiler
recognises. </P>
<P>My <I>make</I> file for Gnu GCC on Unix systems is for use with
<TT>gmake</TT> rather than <TT>make</TT>. I assume your compiler recognises the
.cpp extension. Ordinary <TT>make</TT> works with it on the Sun but did not the
Silicon Graphics or HP machines when I had access to them, many years ago. On Linux use <TT>make</TT>. </P>
<P>My make file for the CC compilers works with the ordinary make. </P>
<P>To compile everything with the CC compiler use </P>
<PRE>   make -f nm_cc.mak
</PRE>

<P>or for the gnu compiler use </P>
<PRE>   make -f nm_gnu.mak
</PRE>

<P>There is a line in the make file for CC <TT>rm -f $*.cxx</TT>. Some systems
won't accept this line and you will need to delete it. In this case, if you
have a bad compile and you are using my scheme for linking .cxx files, you will
need to delete the .cxx file link generated by that compile before you can do
the next one. </P>
<H2><A NAME="custom"></A>2.3 Customising</H2>
<P CLASS="small"><A HREF="#compiler">next</A> - <A HREF="#compiler">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>The file <I>include.h</I> sets a variety of options including several
compiler dependent options. You may need to edit include.h to get the options
you require. If you are using a compiler different from one I have worked with
you may have to set up a new section in include.h appropriate for your
compiler. </P>
<P>Borland, Gnu, Microsoft and Watcom are recognised automatically. If
none of these are recognised a default set of options is used. These are fine
for AT&amp;T, HPUX and Sun C++. If you using a compiler I don't know about you
may have to write a new set of options.</P>
<P>There is an option in include.h for selecting whether you use compiler
supported exceptions, simulated exceptions, or disable exceptions. I now set
<I> compiler supported exceptions</I> as the default. Use the option for
compiler supported exceptions <I>if and only if</I> you have set the option on
your compiler to recognise exceptions. Disabling exceptions sometimes helps
with compilers that are incompatible with my exception simulation scheme. </P>
<P>Activate the appropriate statement to make the element type <i>float</i> or 
<i>double</i>. I suggest you leave it at <i>double</i>.</P>
<P>The option <A HREF="#testing">DO_FREE_CHECK</A> is used for tracking memory
leaks and normally should not be activated. </P>
<P>Activate SETUP_C_SUBSCRIPTS if you want to use traditional C style
<A HREF="#elements">element access</A>. Note that this does <I>not</I> change
the starting point for indices when you are using round brackets for accessing
elements or selecting submatrices. It does enable you to use C style square
brackets. This also activates additional constructors for Matrix, ColumnVector 
and RowVector to simplify use with <i>Numerical Recipes in C++</i>.</P>
<P>Activate <TT>#define use_namespace</TT> if you want to use <A
HREF="#namesp">namespaces</A>. Do this only if you are sure your compiler
supports namespaces. If you do turn this option on, be prepared to turn it off
again if the compiler reports inaccessible variables or the linker reports
missing links. </P>
<P>Activate <TT>#define _STANDARD_</TT> to use the standard names for the
included files and to find the floating point precision data using the floating
point standard. This will work with most recent compilers and is done 
automatically for Borland, Gnu, Intel and Microsoft compilers. </P>
<P>If you haven't defined <TT>_STANDARD_</TT> and are using a compiler that
<I>include.h</I> does not recognise and you want to pick up the floating point
precision data from <I>float.h</I> then activate <TT>#define use_float_h</TT>.
Otherwise the floating point precision data will be accessed from
<I>values.h</I>. You may need to do this with computers from Digital, in
particular.</P>
<P>There is a line</P>
<pre>   //#define set_unix_options</pre>
<p>You can activate this if you are using a Linux or Unix system. It is not used 
by <i>Newmat</i> but is used by some of my other programs.</p>
<H2><A NAME="compiler"></A>2.4 Compilers</H2>
<P CLASS="small"><A HREF="#atandt">next</A> - <A HREF="#update">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH="100%">
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%">
<A HREF="#atandt">2.4.1 AT&amp;T </A><BR>
<A HREF="#borland">2.4.2 Borland </A><BR>
<A HREF="#gcc">2.4.3 Gnu G++ </A><BR>
<A HREF="#hpux">2.4.4 HPUX </A> </TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%">
<A HREF="#intel">2.4.5 Intel </A><BR>
<A HREF="#microso">2.4.6 Microsoft </A><BR>
<A HREF="#sun">2.4.7 Sun </A><BR>
<A HREF="#watcom">2.4.8 Watcom </A></TD>
</TR>
</TABLE>
<P>I have tested this library on a number of compilers. Here are the levels of
success and any special considerations. In most cases I have chosen code that
works under all the compilers I have access to, but I have had to include some
specific work-arounds for some compilers. For the newest PC versions,  I use a Pentium 
4 computer running windows XP or Linux
(Red Hat workstation version). The Unix versions are on a Sun Sparc
station. Thanks to Victoria University for access to the Sparc. </P>
<P>I have set up a block of code for each of the compilers in include.h. Turbo,
Borland, Gnu, Microsoft and Watcom are recognised automatically. There is a
default option that works for AT&amp;T, Sun C++ and HPUX. So you don't
have to make any changes for these compilers. Otherwise you may have to build
your own set of options in include.h. </P>
<H2><A NAME="atandt"></A>2.4.1 AT&amp;T</H2>
<P CLASS="small"><A HREF="#borland">next</A> - <A HREF="#borland">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<P>The AT&amp;T compiler used to be available on a wide variety of Unix
workstations. I don't know if anyone still uses it. However the AT&amp;T options are 
the default if your compiler is not recognised.</P>
<P>AT&amp;T C++ 2.1; 3.0.1 on a Sun: Previous versions worked on these
compilers, which I no longer have access to. </P>
<P>In AT&amp;T 2.1 you may get an error when you use an expression for the
single argument when constructing a Vector or DiagonalMatrix or one of the
Triangular Matrices. You need to evaluate the expression separately. </P>
<H2><A NAME="borland"></A>2.4.2 Borland</H2>
<P CLASS="small"><A HREF="#gcc">next</A> - <A HREF="#gcc">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<H3>Newer compilers</H3>
<P><b>Borland Builder version 8.</b> My tests here have been compiling from a 
make file using <i>nm_b58.mak </i> and making a console program. This works fine 
except that my values of <i>LnMinimum</i> and <i>LnMaximum</i> in <i>precisio.h</i> 
are not correct.</P>
<P><b>Borland Builder version 6:</b> My tests have been on the <i>personal</i> 
version. See the notes for version 5. If you are 
compiling with a make file you can use <i>nm_b56.mak</i> as a model. You can set 
the <i>newmat</i> options to use namespace and the standard library. If you are 
compiling a GUI program you may need to comment out the line defining <i>
TypeDefException</i> in <i>include.h. </i>I don't believe exceptions work 
completely correctly in either version 5 or version 6. However, this does not 
seem to be a problem with my use of them in <i>newmat</i>.</P>
<P><B>Borland Builder version 5:</B> This works fine in console mode and no
special editing of the source codes is required. I haven't tested it in GUI
mode. You can set the <i>newmat</i> options to use namespace and the standard library. <b> You
should turn <I>off</I> the Borland option to use pre-compiled headers.</b> There 
are notes on compiling with the IDE on my <a href="#where">website</a>. 
Alternatively you can use the <i>nm_b55.mak</i> make file<i>.</i></P>
<P><B>Borland Builder version 4</B>: I have successfully used this on older 
versions of newmat using the
console wizard (menu item file/new - select new tab). Use compiler
exceptions. Suppose you are compiling my test program <I>tmt</I>. Rename my
<I>main()</I> function in <I>tmt.cpp</I> to <I>my_main()</I>. Rename
<I>tmt.cpp</I> to <I>tmt_main.cpp</I>. Borland will generate a new file
<I>tmt.cpp</I> containing their <I>main()</I> function. Put the line <TT>int
my_main();</TT> above this function and put <TT>return my_main();</TT> into the
body of <I>main()</I>.</P>
<P><B>Borland compiler version 5.5</B>: this is the free C++ compiler available
from Borland's web site. I suggest you use
the compiler supported exceptions and turn on <I>standard</I> in include.h. You
can use the make file <i>nm_b55.mak</i> after editing to correct the file locations for
your system.</P>
<H3>Older compilers</H3>
<P><B>Borland C++  5.02</B>:</P>
<P>I am not longer checking compatibility with this compiler.</P>
<P>Use the large or 32 bit flat model. If you are not debugging, turn off the 
options that collect debugging information. Use my simulated exceptions.</P>
<P>When running my test program under ms-dos you may run out of memory. Either
compile the test routine to run under <I>easywin</I> or use simulated exceptions 
rather than the built in exceptions. </P>
<P>If you can, upgrade to windows 95 or window NT and use the 32 bit console
model. </P>
<P>If you are using the 16 bit large model, don't forget to keep all matrices
less than 64K bytes in length (90x90 for a rectangular matrix if you are using
<TT>double</TT> as your element type). Otherwise your program will crash
without warning or explanation. You will need to break the <A
HREF="#testing">tmt</A> set of test files into several parts to get the program 
to fit into your computer and run without stack overflow. </P>
<P>You can generate make files for versions 5  with my <a href="genmake.htm">
genmake</a> utility.</P>
<P><b>Borland C++ 3 and 4</b>.</P>
<P>The program used to compile in version 3.1 if you enable the <i>simulated booleans</i> 
- comment <i>out</i> the line <tt>#define bool_LIB 0</tt> in <i>include.h</i> 
and use the <i>simulated exceptions</i>. I haven't checked the latest versions 
of Newmat. The main test program is too large to run 
unless you break it up into several parts. I haven't tried it under version 4.</P>
<H2><A NAME="gcc"></A>2.4.3 Gnu G++</H2>
<P CLASS="small"><A HREF="#hpux">next</A> - <A HREF="#hpux">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<P><B>Gnu G++ 3, 4 (Linux), 3 (Sun):</B> These work OK. If you are using a much earlier version
see if you can upgrade. It&nbsp; used to work with 2.95 and 2.96 but I don't 
have access to these now. You can't use <i>standard</i> with the 2.9X 
versions. The namespace option worked with 2.96 on Linux but not with 2.95 on 
the Sun. Standard is automatically turned on with the 3.X.</P>
<P>This version of Newmat is not compatible with versions 2.6 or earlier.</P>
<H2><A NAME="hpux"></A>2.4.4 HP-UX</H2>
<P CLASS="small"><A HREF="#intel">next</A> - <A HREF="#intel">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<P>HP 9000 series HP-UX. I no longer have access to this compiler. Newmat09
worked without problems with the simulated exceptions; haven't tried the
built-in exceptions. </P>
<P>With recent versions of the compiler you may get warning messages like
<TT>Unsafe cast between pointers/references to incomplete classes</TT>. At
present, I think these can be ignored. </P>

<h2><a name="intel"></a>2.4.5 Intel</h2>
<P CLASS="small"><A HREF="#microso">next</A> - <A HREF="#microso">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>

<P>Newmat works correctly with the Intel 10 C++ compiler for Windows and Linux. (Not tested for the other versions). Standard is 
automatically turned on for both the Linux versions and Windows versions. If 
this causes a problem for the version you are using you can find the lines in <i>
include.h</i> that control this and comment them out. Note that the Intel compiler for 
Linux is <i>free</i> for non-commercial use. (One of the versions of 8.1 gave a 
warning message every time I had something like
<tt>Real x; ... if (x==0.0) ...</tt>, which was often. This is now seems to be fixed.)</P>

<H2><A NAME="microso"></A>2.4.6 Microsoft</H2>
<P CLASS="small"><A HREF="#sun">next</A> - <A HREF="#sun">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<H3>Newer versions</H3>
<P>See my <A HREF="#where">web site</A> for instructions how to work
Microsoft's IDE.<B></B></P>
<P><b>Microsoft Visual C++ 7, 7.1, 8:</b> These work OK. All my tests have 
been in console mode. You can turn on my namespace option. Standard is turned on 
by default for these versions. </P>
<P><B>Microsoft Visual C++ 6</B>: <b>Get the latest service pack</b>. I have tried this 
in console mode and it seems to work satisfactorily. Use the compiler supported exceptions. You may be able to
use the namespace and standard options. If you want to work under MFC
you may need to <TT>#include &quot;stdafx.h&quot;</TT> at the beginning of each .cpp file 
(or turn off precompiled headers). </P>
<P><B>Microsoft Visual C++ 5</B>: I have tried this in console mode on previous 
versions of Newmat. It
seems to work satisfactorily. There may be a problem with <A
HREF="#namesp">namespace</A> (fixed by Service Pack 3?). <B>Turn optimisation
off</B>. Use the compiler supported exceptions. If
you want to work under MFC&nbsp;
you may need to <TT>#include &quot;stdafx.h&quot;</TT> at the
beginning of each .cpp file (or turn off precompiled headers).</P>
<H3>Older versions</H3>
<p>I doubt whether these will work.</p>
<H2><A NAME="sun"></A>2.4.7 Sun</H2>
<P CLASS="small"><A HREF="#watcom">next</A> - <A HREF="#watcom">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<P><B>Sun C++ (version = ?):</B> This seems to work fine with 
compiler supported exceptions. <B>Sun C++ (version
5):</B> There was a problem with exceptions. If you use my simulated
exceptions the non-linear optimisation programs hang. If you use the compiler
supported exceptions my tmt and test_exc programs crash. You should
<I>disable</I> exceptions.</P>
<H2><A NAME="watcom"></A>2.4.8 Watcom</H2>
<P CLASS="small"><A HREF="#update">next</A> - <A HREF="#update">skip</A> -
<A HREF="#compiler">up</A> - <A HREF="#top">start</A></P>
<P><b>Open Watcom (version 1.7a):</b> this works. You can set the standard 
option in <tt>include.h</tt>. The <i>scientific</i> and <i>fixed</i> 
manipulators don't work.</P>
<P><B>Watcom C++ (version 10a):</B> this used to work, I don't know if it works 
now. </P>
<H2><A NAME="update"></A>2.5 Updating from previous
versions</H2>
<P CLASS="small"><a href="#except_1">next</a> - <a href="#except_1">skip</a> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P><b>Newmat11 </b>- if you are upgrading from earlier versions note the 
following:</P>
<ul>
  <li>The class <i>Exception</i> in <i>myexcept.h</i> has been replaced by <i>
  BaseException</i> and a <b>typedef</b> statement included so programs that use 
  the <i>Exception</i> class will still work. If the <i>Exception</i> class was 
  causing a problem comment out the line defining <i>TypeDefException</i> in <i>
  include.h</i>.</li>
  <li>Newmat11 does not support the TEMP_DESTROYED_QUICKLY options so won't work 
	with very old versions of Gnu G++</li>
  <li>Old AT&amp;T work-arounds are removed</li>
  <li>The simulated Booleans class is now stored at the end of include.h. In 
  general, I am not testing with compilers that don't support <b>bool</b>.</li>
  <li>QRZ and QRZT no longer throw an exception if they generate a singular 
	triangular matrix</li>
  <li>Converting functions to lower case. Both the old and new names will work. See the list of <a href="#function">
  functions</a> for new and old names</li>
	<li><i>nm_i5.mak</i>, <i>nm_il5.mak</i> replaced by <i>nm_i8.mak</i>, <i>
	nm_il8.mak</i></li>
	<li>extra file <i>nm_misc.cpp</i> to be included in compilation</li>
	<li><i>scientific</i> and <i>fixed</i> manipulators now recognised in output 
	statements. If you want <i>fixed</i> output and have been using <i>
	scientific</i> in a previous output statement, you may need to include the
	<i>fixed</i> manipulator. (Doesn't work in Visual C++, version 6 and open 
	Watcom).</li>
</ul>
<P><B>Newmat10</B> includes new <A HREF="#scalar2">maxima, minima</A>,
<A HREF="#scalar3">determinant, dot product and Frobenius norm</A> functions, a
faster <A HREF="#fft">FFT</A>, revised <A HREF="#make">make</A> files for GCC
and CC compilers, several corrections, new <A HREF="#dimen">ReSize</A> function, <a href="#constr">IdentityMatrix</a> 
and <a href="#binary">Kronecker Product</a>. Singular values from <a href="#svd">
SVD</a> are sorted. The program files include a new file, <TT>newfft.cpp</TT>, so you
will need to include this in the list of files in your IDE and make files. There 
is also a new test file tmtm.cpp.
<A HREF="#pointer">Pointer arithmetic</A> now mostly meets requirements of
standard. You can use <A HREF="#entering">&lt;&lt;</A> to load data into rows
of a matrix. The <A HREF="#custom">default options</A> in include.h have been
changed. If you are updating from a beta version of newmat09 look through the
next section as there were some late changes to newmat09. </P>
<P>If you are upgrading from <B>newmat08</B> note the following:
</P>
<UL>
<LI>Boolean, TRUE, FALSE are now bool, true, false. See
<A HREF="#custom">customising</A> if your compiler supports the bool class.
</LI>
<LI>ReDimension is now <A HREF="#dimen">ReSize</A>.  
</LI>
<LI>The <A HREF="#except">simulated exception</A> package has
been updated. </LI>
<LI>Operators <TT>==</TT>, <TT>!=</TT>, <TT>+=</TT>,
<TT>-=</TT>, <TT>*=</TT>, <TT>|=</TT>, <TT>&amp;=</TT> are now supported as
<A HREF="#binary">binary</A> matrix operators. </LI>
<LI><TT>A+=f</TT>, <TT>A-=f</TT>, <TT>A*=f</TT>, <TT>A/=f</TT>,
<TT>f+A</TT>, <TT>f-A</TT>, <TT>f*A</TT> are supported for <A HREF="#matscal">A
matrix, f scalar</A>. </LI>
<LI><A HREF="#trigtran">Fast trigonometric transforms</A>.
</LI>
<LI><A HREF="#unary">Reverse</A> function for reversing order of
elements in a vector or matrix. </LI>
<LI><A HREF="#scalar3">IsSingular</A> function. </LI>
<LI>An option is included for defining <A
HREF="#namesp">namespaces</A>. </LI>
<LI>Dummy inequality operators are defined for compatibility
with the STL. </LI>
<LI>The row/column classes in newmat3.cpp have been modified to
improve efficiency and correct an invalid use of pointer arithmetic. Most users
won't be using these classes explicitly; if you are, please contact me for
details of the changes. </LI>
<LI>Matrix LU decomposition rewritten (faster for large arrays).
</LI>
<LI>The sort function rewritten (faster). </LI>
<LI>The documentation files newmata.txt and newmatb.txt have
been amalgamated and both are included in the hypertext version. </LI>
<LI>Some of the <A HREF="#make">make</A> files reorganised
again. </LI>
</UL>
<P>If you are upgrading from <B>newmat07</B> note the following:
</P>
<UL>
<LI>.cxx files are now .cpp files. Some versions of won't accept
.cpp. The <I>make</I> files for Gnu and AT&amp;T link the .cpp files to .cxx
files before compilation and delete the links after compilation.  </LI>
<LI>An <A HREF="#custom">option</A> in include.h allows you to
use compiler supported exceptions, simulated exceptions or disable exceptions.
Edit the file include.h to select one of these three options. Don't simulate
exceptions if you have set your compiler's option to implement exceptions.
</LI>
<LI>New <A HREF="#qr">QR decomposition</A> functions.  
</LI>
<LI>A <A HREF="#nonlin">non-linear least squares</A> class.
</LI>
<LI>No need to explicitly set the AT&amp;T option in include.h.
</LI>
<LI><A HREF="#binary">Concatenation and elementwise
multiplication</A>. </LI>
<LI>A new <A HREF="#unspec">GenericMatrix</A> class.  
</LI>
<LI><A HREF="#scalar3">Sum</A> function. </LI>
<LI>Some of the <A HREF="#make">make</A> files reorganised.
</LI>
</UL>
<P>If you are upgrading from <B>newmat06</B> note the following:
</P>
<UL>
<LI>If you are using &lt;&lt; to load a Real into a submatrix
change this to =. </LI>
</UL>
<P>If you are upgrading from <B>newmat03</B> or <B>newmat04</B>
note the following </P>
<UL>
<LI>.hxx files are now .h files </LI>
<LI>real changed to Real </LI>
<LI>BOOL changed to Boolean </LI>
<LI>CopyToMatrix changed to AsMatrix, etc </LI>
<LI>real(A) changed to A.AsScalar() </LI>
</UL>
<P>The current version is quite a bit longer that newmat04, so
if you are almost out of space with newmat04, don't throw newmat04 away until
you have checked your program will work under this version. </P>
<P>See the <A HREF="#changes">change history</A> for other changes. </P>
<h2><a name="except_1"></a>2.6 Catching exceptions</h2>
<P CLASS="small"><A HREF="#example">next</A> - <A HREF="#example">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>This section applies particularly to people using <i>compiler supported</i>
exceptions rather than my <i>simulated</i> exceptions. </P>
<hr>
<P><b>If newmat detects an error it will throw an exception. It is important that
you catch this exception and print the error message</b>. <b>Otherwise you will get an
unhelpful message like <i>abnormal termination</i>.</b></P>
<hr>
<P>I suggest you set up your
main program like&nbsp; </P>
<pre>#define WANT_STREAM             // or #include &lt;iostream&gt;
#include &quot;newmat.h&quot;             // or #include &quot;newmatap.h&quot;
#include &quot;newmatio.h&quot;           // if you are using my matrix output functions

main()
{
   try
   {
      ... your program here
   }
   // catch exceptions thrown by my programs
   catch(BaseException) { cout &lt;&lt; BaseException::what() &lt;&lt; endl; }
   // catch exceptions thrown by other people's programs
   catch(...) { cout &lt;&lt; &quot;exception caught in main program&quot; &lt;&lt; endl; }
   return 0;
}</pre>
<P>Or see my file <i>nm_ex1.cpp</i> for an easy way of organising this.</P>
<P>If you are using a GUI version rather a console version of the program you 
will need to catch the exception and display the error message in a pop-up 
window.</P>
<P>If you are using my simulated exceptions or have set the disable exceptions
option in <i>include.h</i> then uncaught exceptions automatically print the
error message generated by the exception so you can ignore this section.
Alternatively use <i>Try</i>, <i>Catch</i> and <i>CatchAll</i> in place of <i>try</i>, <i>catch</i>
and <i>catch(...)</i> in the preceding code.</P>
<P>See the <a href="#except">section on exceptions</a> for more information on
the exception structure. </P>
<H2><A NAME="example"></A>2.7 Examples</H2>
<P CLASS="small"><A HREF="#testing">next</A> - <A HREF="#testing">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>I include a number of example files. See the sections on <a href="#make">make</a> 
files and on <a href="#compiler">compilers</a> for information about compiling 
them.</P>
<P><b>Invert matrix:</b> <tt>nm_ex1.cpp</tt>. Load values into a 4x4 matrix; invert it 
and check the result. The output is in <tt>nm_ex1.txt</tt>.</P>
<P><b>Eigenvalues and eigenvectors of Hilbert matrix:</b> <tt>nm_ex2.cpp</tt>. Calculate the 
eigenvalues and eigenvectors of a 7x7 Hilbert matrix. The output is in 
<tt>nm_ex2.txt</tt>.</P>
<P>Values in <tt>precisio.h</tt>: <tt>nm_ex3.cpp</tt>.</P>
<P><b>Linear regression example</b>: <TT>example.cpp</TT>. This gives a  linear
regression example using five different algorithms. The correct output is given
in <TT>example.txt</TT>. The program carries out a rough check that no memory
is left allocated on the heap when it terminates. See the section on
<A HREF="#testing">testing</A> for a comment on the reliability of this check
and the use of the DO_FREE_CHECK option. </P>
<P>Other example files are <TT>nl_ex.cpp</TT> and <TT>garch.cpp</TT> for
demonstrating the non-linear fitting routines, <TT>sl_ex</TT> for demonstrating
the solve function and <TT>test_exc</TT> for demonstrating and testing exception 
handling. </P>
<H2><A NAME="testing"></A>2.8 Testing</H2>
<P CLASS="small"><A HREF="#bugs">next</A> - <A HREF="#bugs">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<P>The library package contains a comprehensive test program in the form of a
series of files with names of the form tmt?.cxx. The files consist of a large
number of matrix formulae all of which evaluate to zero (except the first one
which is used to check that we are detecting non-zero matrices). The printout
should state that it has found just one non-zero matrix. </P>
<P>The test program should be run with <I>Real</I> typedefed to <I>double</I>
rather than <I>float</I> in <A HREF="#custom">include.h</A>. </P>
<P>Make sure the <A HREF="#elements">C subscripts</A> are enabled if you want
to test these. </P>
<P>If you are carrying out some form of bounds checking, for example, with
Borland's <I>CodeGuard</I>, then disable the testing of the <A
HREF="#nric">Numerical Recipes in C</A> interface. Activate the statement
<TT>#define DONT_DO_NRIC</TT> in tmt.h. </P>
<P>Various versions of the make file (extension .mak) are included with the
package. See the section on <A HREF="#make">make files</A>. </P>
<P>The program also allocates and deletes a large block and small block of
memory before it starts the main testing and then at the end of the test. It
then checks that the blocks of memory were allocated in the same place. If not,
then one suspects that there has been a memory leak. i.e. a piece of memory has
been allocated and not deleted.</P>
<P>This is not  foolproof. For example, programs may allocate extra print buffers
while the program is running. I have tried to overcome this by doing a print
before I allocate the first memory block. Programs may allocate memory for
different sized items in different places, or might not allocate items
consecutively. Or they might mix the items with memory blocks from other
programs. Nevertheless, I seem to get consistent answers from <i>some</i> of the
compilers I work with, so I think this is a worthwhile test. The compilers that 
the test seems to work for include the Borland compilers, Microsoft VC++ 6 , 
Watcom, and Gnu 2.96 for Linux.</P>
<P>If the <A HREF="#custom">DO_FREE_CHECK</A> option in include.h is activated,
the program checks that each <TT>new</TT> is balanced with exactly one
<TT>delete</TT>. This provides a more definitive test of no memory leaks. There
are additional statements in myexcept.cpp which can be activated to print out
details of the memory being allocated and released. </P>
<P>I have included a facility for checking that each piece of code in the
library is really exercised by the test routines. Each block of code in the
main part of the library contains a word <TT>REPORT</TT>. <TT>newmat.h</TT> has
a line defining <TT>REPORT</TT> that can be activated (deactivate the dummy
version). This gives a printout of the number of times each of the
<TT>REPORT</TT> statements in the <TT>.cpp</TT> files is accessed. Use a grep
with line numbers to locate the lines on which <TT>REPORT</TT> occurs and
compare these with the lines that the printout shows were actually accessed.
One can then see which lines of code were not accessed. </P>
<H2><A NAME="bugs"></A>2.9 Bugs</H2>
<P CLASS="small"><a href="#problemareas">next</a> - <a href="#problemareas">skip</a> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A> </P>
<UL>
<LI>Small memory leaks may occur when an exception is thrown and caught. </LI>
<LI>My exception scheme may not be not properly linked in with the standard
library exceptions. In particular, my scheme may fail to catch out-of-memory
exceptions.</LI>
</UL>
<h2><a name="problemareas"></a>2.10 Problem areas</h2>
<P CLASS="small"><A HREF="#files">next</A> - <A HREF="#files">skip</A> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A> </P>
<p>This section lists parts of <i>Newmat</i> which users (including me) have found 
difficult or unnatural. Also see the <a href="nmfaq.htm">newmat FAQ list</a> on 
my <a href="#where">website</a>.</p>
<p><b>Invert, element access, matrix multiply etc causes the program to crash</b></p>
<blockquote>
<p CLASS="small">Newmat throws an exception when it detects an error. This can cause a program 
crash unless the exception is caught with a <i>catch</i> statement. See
<a href="#except_1">catching exceptions</a>.</p>
</blockquote>
<p><b><tt>1x1</tt> matrix not automatically converted to a Real</b></p>
<blockquote>
<p CLASS="small">Use the <a href="#ch_type">as_scalar()</a> member function or the
<a href="#scalar3">dotproduct()</a> function to take the dot 
product of two vectors.</p>
</blockquote>
<p><b>Constructors do not initialise elements</b></p>
<blockquote>
<p CLASS="small">For example, <tt>Matrix A(4,5);</tt> does not initialise the elements of <tt>A</tt>. Use the 
statement <tt>A=0.0</tt> to set the values to zero.</p>
</blockquote>
<p><b><i>resize</i> does not initialise elements</b></p>
<blockquote>
<p CLASS="small">For example, <tt>A.resize(5,6);</tt> does not set the elements of <tt>A</tt>. 
If you want to keep values use <tt>resize_keep</tt>. See <a href="#dimen">resize</a>.</p>
</blockquote>
<p><b>Setting Matrix to a scalar sets all the values</b></p>
<blockquote>
<p CLASS="small"><tt>A(1,3) = 0.0;</tt> sets one element of a Matrix to zero. <tt>A = 0.0;</tt> sets all the 
elements to zero. This is very convenient but also a source of error that is 
 
hard to see if you wanted <tt>A(1,3) = 0.0;</tt> but left out the element details.</p>
</blockquote>
<p><b>Symmetry not detected automatically</b></p>
<blockquote>
<p CLASS="small">For example, <tt>SymmetricMatrix SM = A.t() * A;</tt> will fail. Use <tt>SymmetricMatrix 
SM; SM &lt;&lt; A.t() * A;</tt></p>
</blockquote>
<p><b><tt>&lt;&lt;</tt> does not work with constructors</b></p>
<blockquote>
<p CLASS="small">For example, <tt>SymmetricMatrix SM &lt;&lt; A.t() * A;</tt> does not work. </p>
</blockquote>
<p><b>Multiple multiplication may be inefficient</b></p>
<blockquote>
<p CLASS="small">For example, <tt>A * B * X</tt> where <tt>A</tt> and <tt>B</tt> are matrices and <tt>X</tt> is a column vector is 
likely to be much slower than <tt>A * (B * X)</tt>.</p>
</blockquote>
<p>&nbsp;</p>
<H2><A NAME="files"></A>2.11 List of files</H2>
<P CLASS="small"><a href="#refer">next</a> - <a href="#refer">skip</a> -
<A HREF="#starting">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH="100%" BORDER="0" cellspacing="0" style="border-collapse: collapse" bordercolor="#111111" cellpadding="0">
<TR>
<TD width="25%"><B>Documentation</B></TD>
<TD width="25%">readme.txt</TD>
<TD width="50%">readme file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm11.htm</TD>
<TD width="50%">documentation file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">add_time.pgn</TD>
<TD width="50%">image used by nm11.htm</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">rbd.css</TD>
<TD width="50%">style sheet for nm11.htm</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">_newmat.dox</TD>
<TD width="50%">description file for Doxygen</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">_rbd_com.dox</TD>
<TD width="50%">description file for Doxygen</TD>
</TR>
<TR>
<TD width="25%"><b>Header files</b></TD>
<TD width="25%">controlw.h</TD>
<TD width="50%">control word definition file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">include.h</TD>
<TD width="50%">details of include files and options</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">myexcept.h</TD>
<TD width="50%">general exception handler definitions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat.h</TD>
<TD width="50%">main matrix class definition file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatap.h</TD>
<TD width="50%">applications definition file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatio.h</TD>
<TD width="50%">input/output definition file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatnl.h</TD>
<TD width="50%">non-linear optimisation definition file</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatrc.h</TD>
<TD width="50%">row/column functions definition files</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatrm.h</TD>
<TD width="50%">rectangular matrix access definition files</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">precisio.h</TD>
<TD width="50%">numerical precision constants</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">solution.h</TD>
<TD width="50%">one dimensional solve definition file</TD>
</TR>
<TR>
<TD width="25%"><B>Program files</B></TD>
<TD width="25%">bandmat.cpp</TD>
<TD width="50%">band matrix routines</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">cholesky.cpp</TD>
<TD width="50%">Cholesky decomposition</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">evalue.cpp</TD>
<TD width="50%">eigenvalues and eigenvector calculation</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">fft.cpp</TD>
<TD width="50%">fast Fourier, trig. transforms</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">hholder.cpp</TD>
<TD width="50%">QR routines</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">jacobi.cpp</TD>
<TD width="50%">eigenvalues by the Jacobi method</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">myexcept.cpp</TD>
<TD width="50%">general error and exception handler</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newfft.cpp</TD>
<TD width="50%">new fast Fourier transform</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat1.cpp</TD>
<TD width="50%">type manipulation routines</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat2.cpp</TD>
<TD width="50%">row and column manipulation functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat3.cpp</TD>
<TD width="50%">row and column access functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat4.cpp</TD>
<TD width="50%">constructors, resize, utilities</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat5.cpp</TD>
<TD width="50%">transpose, evaluate, matrix functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat6.cpp</TD>
<TD width="50%">operators, element access</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat7.cpp</TD>
<TD width="50%">invert, solve, binary operations</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat8.cpp</TD>
<TD width="50%">LU decomposition, scalar functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat9.cpp</TD>
<TD width="50%">output routines</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatex.cpp</TD>
<TD width="50%">matrix exception handler</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatnl.cpp</TD>
<TD width="50%">non-linear optimisation</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmatrm.cpp</TD>
<TD width="50%">rectangular matrix access functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_misc.cpp</TD>
<TD width="50%">miscellaneous classes, functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">sort.cpp</TD>
<TD width="50%">sorting functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">solution.cpp</TD>
<TD width="50%">one dimensional solve</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">submat.cpp</TD>
<TD width="50%">submatrix functions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">svd.cpp</TD>
<TD width="50%">singular value decomposition</TD>
</TR>
<TR>
<TD width="25%"><B>Example files</B></TD>
<TD width="25%">nm_ex1.cpp</TD>
<TD width="50%">simple example - invert matrix</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ex1.txt</TD>
<TD width="50%">output from nm_ex1.cpp</TD>
</TR>
<tr>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ex2.cpp</TD>
<TD width="50%">simple example - eigenvalues of Hilbert matrix</TD>
</tr>
<tr>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ex2.txt</TD>
<TD width="50%">output from nm_ex2.cpp</TD>
</tr>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ex3.cpp</TD>
<TD width="50%">scientific format and constants from precisio.h</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ex3.txt</TD>
<TD width="50%">output from nm_ex3.cpp</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">example.cpp</TD>
<TD width="50%">example of use of package</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">example.txt</TD>
<TD width="50%">output from example</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">sl_ex.cpp</TD>
<TD width="50%">example of OneDimSolve routine</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">sl_ex.txt</TD>
<TD width="50%">output from example</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nl_ex.cpp</TD>
<TD width="50%">example of non-linear least squares</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nl_ex.txt</TD>
<TD width="50%">output from example</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">garch.cpp</TD>
<TD width="50%">example of maximum-likelihood fit</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">garch.dat</TD>
<TD width="50%">data file for garch.cpp</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">garch.txt</TD>
<TD width="50%">output from example</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">test_exc.cpp</TD>
<TD width="50%">demonstration exceptions</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">test_exc.txt</TD>
<TD width="50%">output from test_exc.cpp</TD>
</TR>
<TR>
<TD width="25%"><B>Test files</B></TD>
<TD width="25%">tmt.h</TD>
<TD width="50%">header file for test files</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">tmt*.cpp</TD>
<TD width="50%">test files (see the section on <A
HREF="#testing">testing</A>)</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">tmt.txt</TD>
<TD width="50%">output from test files</TD>
</TR>
<TR>
<TD width="25%"><B>Make files</B></TD>
<TD width="25%">nm_gnu.mak</TD>
<TD width="50%">make file for Gnu G++</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_cc.mak</TD>
<TD width="50%">make file for AT&amp;T, Sun and HPUX</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_b55.mak</TD>
<TD width="50%">make file for Borland C++ 5.5</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_b56.mak</TD>
<TD width="50%">make file for Borland Builder C++ 6</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_b58.mak</TD>
<TD width="50%">make file for Borland Builder C++ 8</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_m6.mak</TD>
<TD width="50%">make file for VC++ 6&amp;7</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_m8.mak</TD>
<TD width="50%">make file for VC++ 8</TD>
</TR>
<tr>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_i8.mak</TD>
<TD width="50%">make file for Intel C++ 8,9 under Windows</TD>
</tr>
<tr>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_i10.mak</TD>
<TD width="50%">make file for Intel C++ 10 under Windows</TD>
</tr>
<tr>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_il8.mak</TD>
<TD width="50%">make file for Intel C++ 8,9,10 under Linux</TD>
</tr>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_ow.mak</TD>
<TD width="50%">make file for Open Watcom</TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">newmat.lfl</TD>
<TD width="50%">library file list for use with <a href="genmake.htm">genmake</a></TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">nm_targ.txt</TD>
<TD width="50%">target file list for use with <a href="genmake.htm">genmake</a></TD>
</TR>
<TR>
<TD width="25%">&nbsp;</TD>
<TD width="25%">Makefile.in</TD>
<TD width="50%">used for compiling with <a href="opt.htm">Opt++</a></TD>
</TR>
</TABLE>

<p>&nbsp;</p>
<p>These are the .cpp files you need to include in your make file or project:</p>
<H4>Basic newmat</H4>
<TABLE BORDER="0" CELLPADDING="5" id="table2">
<TR>
<TD VALIGN="top"><UL>
<LI>bandmat.cpp</LI>
<LI>myexcept.cpp</LI>
<LI>newmat1.cpp</LI>
<LI>newmat2.cpp</LI>
</UL>
</TD>
<TD VALIGN="top"><UL>
<LI>newmat3.cpp</LI>
<LI>newmat4.cpp</LI>
<LI>newmat5.cpp</LI>
<LI>newmat6.cpp</LI>
</UL>
</TD>
<TD VALIGN="top"><UL>
<LI>newmat7.cpp</LI>
<LI>newmat8.cpp</LI>
<LI>newmat9.cpp</LI>
<LI>newmatex.cpp</LI>
</UL>
</TD>
<TD VALIGN="top"><UL>
<LI>newmatrm.cpp</LI>
<LI>submat.cpp</LI>
</UL>
</TD>
</TR>
</TABLE>
<H4>Factorisations etc</H4>
<TABLE BORDER="0" id="table3">
<TR>
<TD VALIGN="top"><UL>
<LI>cholesky.cpp</LI>
<LI>evalue.cpp</LI>
<LI>fft.cpp</LI>
</UL>
</TD>
<TD VALIGN="top">
<ul>
	<li>hholder.cpp</li>
	<li>jacobi.cpp</li>
</ul>
</TD>
<TD VALIGN="top">
<ul>
	<li>sort.cpp</li>
<LI>svd.cpp</LI>
</ul>
<p>&nbsp;</TD>
<TD VALIGN="top"><UL>
<LI>newfft.cpp</LI>
<LI>nm_misc.cpp</LI>
</UL>
</TD>
</TR>
</TABLE>
<H4>Nonlinear routines</H4>
<TABLE BORDER="0" id="table4">
<TR>
<TD><UL>
<LI>newmatnl.cpp</LI>
</UL>
</TD>
<TD><UL>
<LI>solution.cpp</LI>
</UL>
</TD>
<TD></TD>
<TD></TD>
</TR>
</TABLE>
<p>&nbsp;</p>

<H2><A NAME="refer"></A>3. Reference manual</H2>
<P CLASS="small"><A HREF="#constr">next</A> - <A HREF="#error">skip</A> -
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH="100%">
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#constr">3.1
Constructors </A><BR>
<A HREF="#elements">3.2 Accessing elements </A><BR>
<A HREF="#copy">3.3 Assignment and copying </A><BR>
<A HREF="#entering">3.4 Entering values </A><BR>
<A HREF="#unary">3.5 Unary operations </A><BR>
<A HREF="#binary">3.6 Binary operations </A><BR>
<A HREF="#matscal">3.7 Matrix and scalar ops </A><BR>
<A HREF="#scalar1">3.8 Scalar functions - size &amp; shape </A><BR>
<A HREF="#scalar2">3.9 Scalar functions - maximum &amp; minimum </A><BR>
<A HREF="#scalar3">3.10 Scalar functions - numerical </A><BR>
<A HREF="#submat">3.11 Submatrices </A><BR>
<A HREF="#dimen">3.12 Change dimension </A><BR>
<A HREF="#ch_type">3.13 Change type </A><BR>
<A HREF="#solve">3.14 Multiple matrix solve </A><BR>
<A HREF="#memory">3.15 Memory management </A><BR>
<A HREF="#efficien">3.16 Efficiency</A><br>
<A
HREF="#output">3.17 Output</A><br>
<A HREF="#unspec">3.18 Accessing unspecified type</A></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%">  
<A HREF="#cholesky">3.19 Cholesky decomposition </A><BR>
<A HREF="#qr">3.20 QR decomposition </A><BR>
<A HREF="#svd">3.21 Singular value decomposition </A><BR>
<A HREF="#evalues">3.22 Eigenvalue decomposition </A><BR>
<A HREF="#sorting">3.23 Sorting </A><BR>
<A HREF="#fft">3.24 Fast Fourier transform </A><BR>
<A HREF="#trigtran">3.25 Fast trigonometric transforms </A><BR>
<A HREF="#nric">3.26 Numerical recipes in C </A><BR>
<A HREF="#except">3.27 Exceptions </A><BR>
<A HREF="#cleanup">3.28 Cleanup following exception </A><BR>
<A HREF="#nonlin">3.29 Non-linear applications </A><BR>
<A HREF="#stl">3.30 Standard template library </A><BR>
<A HREF="#namesp">3.31 Namespace </A><br>
<a href="#upd_chol">3.32 Updating the Cholesky matrix</a><br>
<a href="#RealStarStar">3.33 Accessing C functions</a><br>
<a href="#SimpleIntArray">3.34 Simple integer array class</a><br>
<a href="#extend">3.35 Extend orthonormal set of columns</a><br>
<a href="#misc_fn">3.36 Miscellaneous functions</a></TD>
</TR>
</TABLE>
<p>See also <a href="#function">function summary</a>.</p>

<H2><A NAME="constr"></A>3.1 Constructors</H2>
<P CLASS="small"><A HREF="#elements">next</A> - <A HREF="#elements">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>To construct an <I>m</I> x <I>n</I> matrix, <TT>A</TT>, (<I>m</I> and
<I>n</I> are integers) use </P>
<PRE>    Matrix A(m,n);
</PRE>

<P>The SquareMatrix, UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and
DiagonalMatrix types are square. To construct an <I>n</I> x <I>n</I> matrix
use, for example </P>
<PRE>    SquareMatrix SQ(n);
    UpperTriangularMatrix UT(n);
    LowerTriangularMatrix LT(n);
    SymmetricMatrix S(n);
    DiagonalMatrix D(n);
</PRE>

<P>Band matrices need to include bandwidth information in their constructors. 
</P>
<PRE>    BandMatrix BM(n, lower, upper);
    UpperBandMatrix UB(n, upper);
    LowerBandMatrix LB(n, lower);
    SymmetricBandMatrix SB(n, lower);
</PRE>

<P>The integers <I>upper</I> and <I>lower</I> are the number of non-zero
diagonals above and below the diagonal (<I>excluding</I> the diagonal)
respectively.&nbsp; The UpperBandMatrix and LowerBandMatrix are upper and lower 
triangular band matrices. So an UpperBandMatrix is essentially a BandMatrix with
<i>lower</i> = 0 and a LowerBandMatrix is a BandMatrix with <i>upper</i> = 0.</P>
<P>The RowVector and ColumnVector types take just one argument in their
constructors: </P>
<PRE>    RowVector RV(n);
    ColumnVector CV(n);
</PRE>

<P><b>These constructors do <EM>not</EM> initialise the elements of the matrices.
</b>
To set all the elements to zero use, for example, </P>
<PRE>    Matrix A(m, n); A = 0.0;
</PRE>

<P>The IdentityMatrix takes one argument in its constructor specifying its 
dimension.</P>
<pre>    IdentityMatrix I(n);</pre>
<p>The value of the diagonal elements <b>is</b> set to 1 by default, but you can 
change this value as with other matrix types. </p>

<P>You can also construct vectors and matrices without specifying the
dimension. For example </P>
<PRE>    Matrix A;
</PRE>

<P>In this case the dimension must be set by an <A HREF="#copy">assignment
statement</A> or a <A HREF="#dimen">resize statement</A>. </P>
<P>You can also use a constructor to set a matrix equal to another matrix or
matrix expression. </P>
<PRE>    Matrix A = UT;
    Matrix A = UT * LT;
</PRE>

<P>Only conversions that don't lose information are supported - eg you cannot
convert an upper triangular matrix into a diagonal matrix using =. </P>
<H2><A NAME="elements"></A>3.2 Accessing elements 
</H2>
<P CLASS="small"><A HREF="#copy">next</A> - <A HREF="#copy">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>Elements are accessed by expressions of the form <TT>A(i,j)</TT> where
<I>i</I> and <I>j</I> run from 1 to the appropriate dimension. Access elements
of vectors with just one argument. Diagonal matrices can accept one or two
subscripts. </P>
<P>This is different from the earliest version of the package in which the
subscripts ran from 0 to one less than the appropriate dimension. Use
<TT>A.element(i,j)</TT> if you want this earlier convention. </P>
<P><TT>A(i,j)</TT> and <TT>A.element(i,j)</TT> can appear on either side of an
= sign. </P>
<P>If you activate the <TT>#define SETUP_C_SUBSCRIPTS</TT> in
<TT>include.h</TT> you can also access elements using the traditional C style
notation. That is <TT>A[i][j]</TT> for matrices (except diagonal) and
<TT>V[i]</TT> for vectors and diagonal matrices. The subscripts start at zero
(i.e. like element) and there is <I>no</I> range checking. Because of the
possibility of confusing <TT>V(i)</TT> and <TT>V[i]</TT>, I suggest you do
<I>not</I> activate this option unless you really want to use it.</P>
<P>Symmetric matrices are stored as lower triangular matrices. It is important
to remember this if you are using the <TT>A[i][j]</TT> method of accessing
elements. Make sure the first subscript is greater than or equal to the second
subscript. However, if you are using the <TT>A(i,j)</TT> method the program
will swap <TT>i</TT> and <TT>j</TT> if necessary; so it doesn't matter if you
think of the storage as being in the upper triangle (but it <I>does</I> matter
in some other situations such as when <A HREF="#entering">entering</A> data).</P>
<P>The IdentityMatrix type does not support element access.</P>
<P>For interfacing with traditional C functions that involve one and two 
dimensional arrays see <a href="#RealStarStar">accessing C functions</a>.</P>
<H2><A NAME="copy"></A>3.3 Assignment and copying</H2>
<P CLASS="small"><A HREF="#entering">next</A> - <A HREF="#entering">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The operator <TT>=</TT> is used for copying matrices, converting matrices,
or evaluating expressions. For example </P>
<PRE>    A = B;  A = L;  A = L * U;
</PRE>

<P>Only conversions that don't lose information are supported. The dimensions
of the matrix on the left hand side are adjusted to those of the matrix or
expression on the right hand side. Elements on the right hand side which are
not present on the left hand side are set to zero. </P>
<P>The operator <TT>&lt;&lt;</TT> can be used in place of <TT>=</TT> where it
is permissible for information to be lost. </P>
<P>For example </P>
<PRE>    SymmetricMatrix S; Matrix A;
    ......
    S &lt;&lt; A.t() * A;
</PRE>

<P>is acceptable whereas </P>
<PRE>    S = A.t() * A;                            // error
</PRE>

<P>will cause a runtime error since the package does not (yet?) recognise
<TT>A.t()*A</TT> as symmetric. </P>
<P>Note that you can <I>not</I> use <TT>&lt;&lt;</TT> with constructors. For
example </P>
<PRE>    SymmetricMatrix S &lt;&lt; A.t() * A;           // error
</PRE>

<P>does <I>not</I> work. </P>
<P>Also note that <TT>&lt;&lt;</TT> cannot be used to load values from a full
matrix into a band matrix, since it will be unable to determine the bandwidth
of the band matrix. </P>
<P>A third copy routine is used in a similar role to <TT>=</TT>. Use </P>
<PRE>    A.inject(D);
</PRE>

<P>to copy the elements of <TT>D</TT> to the corresponding elements of
<TT>A</TT> but leave the elements of <TT>A</TT> unchanged if there is no
corresponding element of <TT>D</TT> (the <TT>=</TT> operator would set them to
0). This is useful, for example, for setting the diagonal elements of a matrix
without disturbing the rest of the matrix. Unlike <TT>=</TT> and
<TT>&lt;&lt;</TT>, inject does not reset the dimensions of <TT>A</TT>, which
must match those of <TT>D</TT>. Inject does <I>not</I> test for no loss of
information. The name <i>Inject</i> can be used instead on <i>inject</i>.</P>
<P>You cannot replace <TT>D</TT> by a matrix expression. The effect of <TT>inject(D)</TT> depends on the type of <TT>D</TT>. If <TT>D</TT> is an
expression it might not be obvious to the user what type it would have. So I
thought it best to disallow expressions. </P>
<P>Inject can be used for loading values from a regular matrix into a band
matrix. (Don't forget to zero any elements of the left hand side that will not
be set by the loading operation).</P>
<P>Both <TT>&lt;&lt;</TT> and inject can be used with submatrix expressions on
the left hand side. See the section on <A HREF="#submat">submatrices</A>. </P>
<P>To set the elements of a matrix to a scalar use operator <TT>=</TT> </P>
<PRE>    Real r; int m,n;
    ......
    Matrix A(m,n); A = r;</PRE>


To swap the values in two matrices <tt>A</tt> and <tt>B</tt> use one of the 
following expressions<pre>   A.swap(B);</pre>
<pre>   swap(A,B);</pre>
The matrices <tt>A</tt> and <tt>B</tt> must be of the same type. This can be any 
of the matrix types listed in the <a href="#constr">section on constructors</a>,
<a href="#solve">CroutMatrix, BandLUMatrix</a> or <a href="#unspec">
GenericMatrix</a>. Swap works by switching pointers and does not do any actual copying 
of the main data arrays.<p>Notes:</p>

<ul>
  <li>When you do a matrix assignment to another matrix or matrix expression with 
either <tt>=</tt> or <tt>&lt;&lt;</tt> the original data array associated with the matrix 
  being assigned to is destroyed 
even if there is no change in length. See the section on <a href="#stor">storage</a>. 
This means, in particular, that pointers to matrix elements - e.g.
<tt>Real* a; a = &amp;(A(1,1));</tt> become invalid. If you want avoid this you 
  can use 
<tt>Inject</tt> rather than <tt>=</tt>. But remember that you may need to zero 
  the matrix first.<br>
  </li>
</ul>

<H2><A NAME="entering"></A>3.4 Entering values</H2>
<P CLASS="small"><A HREF="#unary">next</A> - <A HREF="#unary">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>You can load the elements of a matrix from an array: </P>
<PRE>    Matrix A(3,2);
    Real a[] = { 11,12,21,22,31,33 };
    A &lt;&lt; a;</PRE>
<p>or</p>
<pre>    Matrix A(3,2);
    int a[] = { 11,12,21,22,31,33 };
    A &lt;&lt; a;</pre>
<p>This construction does <I>not</I> check that the numbers of elements match
correctly. This version of <TT>&lt;&lt;</TT> can be used with submatrices on
the left hand side. It is not defined for band matrices.</p>
<p>Note that you enter only the values stored in a matrix. For example</p>
<pre>    SymmetricMatrix A(2);
    Real a[] = { 11,12,22 };
    A &lt;&lt; a;</pre>
<P>Alternatively you can enter short lists using a sequence of numbers
separated by <TT>&lt;&lt;</TT> . </P>
<PRE>    Matrix A(3,2);
    A &lt;&lt; 11 &lt;&lt; 12
      &lt;&lt; 21 &lt;&lt; 22
      &lt;&lt; 31 &lt;&lt; 32;
</PRE>

<P>This does check for the correct total number of entries, although the
message for there being insufficient numbers in the list may be delayed until
the end of the block or the next use of this construction. This does <I>not</I>
work for band matrices or for long lists. It does work for submatrices if the
submatrix is a single complete row. For example </P>
<PRE>    Matrix A(3,2);
    A.row(1) &lt;&lt; 11 &lt;&lt; 12;
    A.row(2) &lt;&lt; 21 &lt;&lt; 22;
    A.row(3) &lt;&lt; 31 &lt;&lt; 32;
</PRE>

<P>Load only values that are actually stored in the matrix. For example </P>
<PRE>    SymmetricMatrix S(2);
    S.row(1) &lt;&lt; 11;
    S.row(2) &lt;&lt; 21 &lt;&lt; 22;
</PRE>

<P>Try to restrict this way of loading data to numbers. You can include
expressions, but these must not call a function which includes the same
construction. </P>

<p>Remember that matrices are stored by rows and that symmetric matrices are
stored as <I> lower</I> triangular matrices when using these methods to enter
data. </p>
<H2><A NAME="unary"></A>3.5 Unary operators</H2>
<P CLASS="small"><A HREF="#binary">next</A> - <A HREF="#binary">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The package supports unary operations </P>
<PRE>    X = -A;           // change sign of elements
    X = A.t();        // transpose
    X = A.i();        // inverse (of square matrix A)
    X = A.reverse();  // reverse order of elements of vector
                      // or matrix (not band matrix)
    ColumnVector X = A.sum_rows();         // sum of elements
                                           // of each row
    RowVector X = A.sum_columns();         // sum of elements
                                           // of each column
    ColumnVector X = A.sum_square_rows();  // sum of squares of
                                           // elements of each row
    RowVector X = A.sum_square_columns();  // sum of squares of
                                           // elements of each column</PRE>

<p class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function name.</p>


<H2><A NAME="binary"></A>3.6 Binary operators</H2>
<P CLASS="small"><A HREF="#matscal">next</A> - <A HREF="#matscal">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The package supports binary operations </P>
<PRE>    X = A + B;       // matrix addition
    X = A - B;       // matrix subtraction
    X = A * B;       // matrix multiplication
    X = A.i() * B;   // equation solve (square matrix A)
    X = A | B;       // concatenate horizontally (concatenate the rows)
    X = A &amp; B;       // concatenate vertically (concatenate the columns)
    X = SP(A, B);    // elementwise product of A and B (Schur product)
    X = KP(A, B);    // Kronecker product of A and B
    X = crossproduct(A, B);          // vector cross product - see notes
    X = crossproduct_rows(A, B);     // cross product of rows
    X = crossproduct_columns(A, B);  // cross product of columns
    bool b = A == B; // test whether A and B are equal
    bool b = A != B; // ! (A == B)
    A += B;          // A = A + B;
    A -= B;          // A = A - B;
    A *= B;          // A = A * B;
    A |= B;          // A = A | B;
    A &amp;= B;          // A = A &amp; B;
    A.SP_eq(B);      // A = SP(A, B);
    &lt;, &gt;, &lt;=, &gt;=     // included for compatibility with STL - see notes
</PRE>

<P>Notes: </P>
<UL>
<LI>If you are doing repeated multiplication. For example <TT>A*B*C</TT>, use
brackets to force the order of evaluation to minimise the number of operations.
If <TT>C</TT> is a column vector and <TT>A</TT> is not a vector, then it will
usually reduce the number of operations to use <TT>A*(B*C)</TT>. </LI>
<LI>In the equation solve example case the inverse is not explicitly
calculated. An LU decomposition of <TT>A</TT> is performed and this is applied
to <TT>B</TT>. This is more efficient than calculating the inverse and then
multiplying. See also <A HREF="#solve">multiple matrix solving</A>. </LI>
<LI>The package does not (yet?) recognise <TT>B*A.i()</TT> as an equation solve
and the inverse of <TT>A</TT> would be calculated. It is probably better to use
<TT>(A.t().i()*B.t()).t()</TT>. </LI>
<LI>Horizontal or vertical concatenation returns a result of type Matrix,
RowVector or ColumnVector. </LI>
<LI>If <TT>A</TT> is <I> m</I> x <I>p</I>, <TT>B</TT> is <I> m</I> x <I>q</I>,
then <TT>A | B</TT> is <I> m</I> x (<I>p</I>+<I>q</I>) with the <I>k</I>-th row
being the elements of the <I>k</I>-th row of <TT>A</TT> followed by the
elements of the <I>k</I>-th row of <TT>B</TT>. </LI>
<LI>If <TT>A</TT> is <I> p</I> x <I>n</I>, <TT>B</TT> is <I> q</I> x <I>n</I>,
then <TT>A &amp; B</TT> is (<I>p</I>+<I>q</I>) x <I> n</I> with the <I>k</I>-th
column being the elements of the <I>k</I>-th column of <TT>A</TT> followed by
the elements of the <I>k</I>-th column of <TT>B</TT>. </LI>
<LI>For complicated concatenations of matrices, consider instead using
<A HREF="#submat">submatrices</A>. </LI>
<LI>See the section on <A HREF="#submat">submatrices</A> on using a submatrix
on the RHS of an expression.</LI>
<LI>crossproduct - assumes A and B are both RowVectors of length 3 or both 
ColumnVectors of length 3. Result is a Matrix of same dimension as A or B (will 
automatically convert to RowVector if A and B are RowVectors and ColumnVector if 
they are ColumnVectors).</LI>
<LI>crossproduct_rows - assumes A and B are of type Matrix with the same number 
of rows and 3 columns. Does a cross product on corresponding pairs of rows.</LI>
<LI>crossproduct_columns - assumes A and B are of type Matrix with&nbsp; the same 
number of columns and 3 rows. Does a cross product on corresponding pairs of 
columns.</LI>
<LI>Two matrices are equal if their difference is zero. They may be of
different types. For the CroutMatrix or BandLUMatrix they must be of the same
type and have all their elements equal. This is not a very useful operator and
is included for compatibility with some container templates. </LI>
<LI>The inequality operators are included for compatibility with some versions 
of the
<A HREF="#stl">standard template library</A>. If actually called, they will
throw an exception. So don't try to sort a <I>list</I> of matrices. </LI>
<LI> A row vector multiplied by a column vector yields a 1x1 matrix, <I>not</I>
a Real. To get a Real result use either <a href="#scalar3">as_scalar or
dotproduct</a>. </LI>
<LI> The result from Kronecker product, <tt>KP(A, B)</tt>, possesses an attribute such as 
upper triangular, lower triangular, band, symmetric, diagonal if both of the 
matrices <TT>A</TT> and <TT>B</TT> have the attribute. If <TT>A</TT> is band and <TT>B</TT> is a square type 
(eg SquareMatrix, Diagonal, LowerTriangularMatrix etc) then 
the result is band.</LI>
<LI> Elementwise product is also known as the Schur product or the Hadamard 
product.</LI>
<LI> See the <a href="#function">function summary list</a> for the older 
depreciated function names.</LI>
</UL>
<P CLASS="small">Remember that the product of symmetric matrices is not
necessarily symmetric so the following code will not run:</P>
<PRE>   SymmetricMatrix A, B;
   .... put values in A, B ....
   SymmetricMatrix C = A * B;   // run time error</PRE>

<P CLASS="small">Use instead</P>
<PRE>   Matrix C = A * B;</PRE>

<P CLASS="small">or, if you <EM>know</EM> the product will be symmetric,
use</P>
<PRE>   SymmetricMatrix C; C &lt;&lt; A * B;
</PRE>

<H2><A NAME="matscal"></A>3.7 Matrix and scalar</H2>
<P CLASS="small"><A HREF="#scalar1">next</A> - <A HREF="#scalar1">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The following expressions multiply the elements of a matrix <TT>A</TT> by a
scalar <TT>f</TT>: <TT>A * f</TT> or <TT>f * A</TT> . Likewise one can divide
the elements of a matrix <TT>A</TT> by a scalar <TT>f</TT>: <TT>A / f</TT> . 
</P>
<P>The expressions <TT>A + f</TT> and <TT>A - f</TT> add or subtract a
rectangular matrix of the same dimension as <TT>A</TT> with elements equal to
<TT>f</TT> to or from the matrix <TT>A</TT> . </P>
<P>The expression <TT>f + A</TT> is an alternative to <TT>A + f</TT>. The
expression <TT>f - A</TT> subtracts matrix <TT>A</TT> from a rectangular matrix
of the same dimension as <TT>A</TT> and with elements equal to <TT>f</TT> . 
</P>
<P>The expression <TT>A += f</TT> replaces <TT>A</TT> by <TT>A + f</TT>.
Operators <TT>-=</TT>, <TT>*=</TT>, <TT>/=</TT> are similarly defined. </P>
<H2><A NAME="scalar1"></A>3.8 Scalar functions of a
matrix - size &amp; shape</H2>
<P CLASS="small"><A HREF="#scalar2">next</A> - <A HREF="#scalar2">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>This page describes functions returning the values associated with the size
and shape of matrices. The following pages describe other scalar matrix
functions. </P>
<PRE>    int m = A.nrows();                     // number of rows
    int n = A.ncols();                     // number of columns
    MatrixType mt = A.type();              // type of matrix
    Real* s = A.data();                    // pointer to array of elements
    const Real* s = A.data();              // pointer to array of elements
                                           //    where A is const
    const Real* s = A.const_data();        // pointer to array of elements
    int l = A.size();                      // length of array of elements
    MatrixBandWidth mbw = A.bandwidth();   // upper and lower bandwidths
</PRE>

<P><TT>MatrixType mt = A.type()</TT> returns the type of a matrix. Use
<TT>mt.value()</TT> to get a string (UT, LT, Rect, Sym, Diag, Band, UB, LB,
Crout, BndLU) showing the type (Vector types are returned as Rect). </P>
<P><TT>MatrixBandWidth</TT> has member functions <TT>upper()</TT> and
<TT>lower()</TT> for finding the upper and lower bandwidths (number of diagonals 
above and below the diagonal, both zero for a diagonal matrix). For non-band 
matrices -1 is returned for both these values.</P>
<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>
<H2><A NAME="scalar2"></A>3.9 Scalar functions of a
matrix - maximum &amp; minimum</H2>
<P CLASS="small"><A HREF="#scalar3">next</A> - <A HREF="#scalar3">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>This page describes functions for finding the maximum and minimum elements
of a matrix. </P>
<PRE>    int i, j;
    Real mv = A.maximum_absolute_value();    // maximum of absolute values
    Real mv = A.minimum_absolute_value();    // minimum of absolute values
    Real mv = A.maximum();                   // maximum value
    Real mv = A.minimum();                   // minimum value
    Real mv = A.maximum_absolute_value1(i);  // maximum of absolute values
    Real mv = A.minimum_absolute_value1(i);  // minimum of absolute values
    Real mv = A.maximum1(i);                 // maximum value
    Real mv = A.minimum1(i);                 // minimum value
    Real mv = A.maximum_absolute_value2(i,j);// maximum of absolute values
    Real mv = A.minimum_absolute_value2(i,j);// minimum of absolute values
    Real mv = A.maximum2(i,j);               // maximum value
    Real mv = A.minimum2(i,j);               // minimum value
</PRE>

<P>All these functions throw an exception if <TT>A</TT> has no rows or no
columns. </P>
<P>The versions <TT>A.maximum_absolute_value1(i)</TT>, etc return the location of
the extreme element in a RowVector, ColumnVector or DiagonalMatrix. The
versions <TT>A.maximum_absolute_value2(i,j)</TT>, etc return the row and column
numbers of the extreme element. If the extreme value occurs more than once the
location of the last one is given. </P>
<P>The versions maximum_absolute_value(A), minimum_absolute_value(A), maximum(A),
minimum(A) can be used in place of A.maximum_absolute_value(),
A.minimum_absolute_value(), A.maximum(), A.minimum(). </P>
<H2><A NAME="scalar3"></A>3.10 Scalar functions of a
matrix - numerical</H2>
<P CLASS="small"><A HREF="#submat">next</A> - <A HREF="#submat">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<PRE>    Real r = A.as_scalar();                // value of 1x1 matrix
    Real ssq = A.sum_square();             // sum of squares of elements
    Real sav = A.sum_absolute_value();     // sum of absolute values
    Real s = A.sum();                      // sum of values
    Real norm = A.norm1();                 // maximum of sum of absolute
                                              values of elements of a column
    Real norm = A.norm_infinity();         // maximum of sum of absolute
                                              values of elements of a row
    Real norm = A.norm_Frobenius();        // square root of sum of squares
                                           // of the elements
    Real t = A.trace();                    // trace
    Real d = A.determinant();              // determinant
    LogAndSign ld = A.log_determinant();   // natural log of determinant
    bool z = A.is_zero();                  // test all elements zero
    bool s = A.is_singular();              // A is a CroutMatrix or
                                              BandLUMatrix
    Real s = dotproduct(A, B);             // dot product of A and B
                                           // interpreted as vectors
</PRE>

<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>

<P><TT>A.log_determinant()</TT> returns a value of type <i>LogAndSign</i>. If ld is of
type <i>LogAndSign</i> use </P>
<PRE>    ld.value()     to get the value of the determinant
    ld.sign()      to get the sign of the determinant (values 1, 0, -1)
    ld.log_value() to get the log of the absolute value.
</PRE>

<P>Note that the direct use of the function <TT>determinant()</TT> will often
cause a floating point overflow exception. </P>
<P><TT>A.is_zero()</TT> returns Boolean value <TT>true</TT> if the matrix
<TT>A</TT> has all elements equal to 0.0. </P>
<P><TT>is_singular()</TT> is defined only for CroutMatrix and BandLUMatrix. It
returns <TT>true</TT> if one of the diagonal elements of the LU decomposition
is exactly zero. </P>
<P><TT>dotproduct(const Matrix&amp; A, const Matrix&amp; B)</TT> converts both
of the arguments to rectangular matrices, checks that they have the same number
of elements and then calculates the first element of <TT>A * </TT>first element
of <TT>B + </TT>second element of <TT>A * </TT>second element of <TT>B +
...</TT> ignoring the row/column structure of <TT>A</TT> and <TT>B</TT>. It is
primarily intended for the situation where <TT>A</TT> and <TT>B</TT> are row or
column vectors. </P>
<P>The versions sum(A), sum_square(A), sum_absolute_value(A), trace(A),
log_determinant(A), determinant(A), norm1(A), norm_infinity(A), norm_Frobenius(A)
can be used in place of A.sum(), A.sum_square(), A.sum_absolute_value(),
A.trace(), A.log_determinant(), A.determinant(A), A.norm1(), A.norm_infinity(), A.norm_Frobenius(). 
</P>
<H2><A NAME="submat"></A>3.11 Submatrices</H2>
<P CLASS="small"><A HREF="#dimen">next</A> - <A HREF="#dimen">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<PRE>    A.submatrix(fr,lr,fc,lc)
</PRE>

<P>This selects a submatrix from <TT>A</TT>. The arguments
<TT>fr</TT>,<TT>lr</TT>,<TT>fc</TT>,<TT>lc</TT> are the first row, last row,
first column, last column of the submatrix with the numbering beginning at 1. 
</P>
<P>I allow <TT>lr = fr-1</TT> or <TT>lc = fc-1</TT> or to indicate that a
matrix of zero rows or columns is to be returned. </P>
<P>A submatrix command may be used in any matrix expression or on the left hand
side of <TT>=</TT>, <TT>&lt;&lt;</TT> or <i>inject</i>. <i>inject</i> does <I>not</I> check
no information loss. You can also use the construction </P>
<PRE>    Real c; .... A.submatrix(fr,lr,fc,lc) = c;
</PRE>

<P>to set a submatrix equal to a constant. </P>
<P>The following are variants of submatrix: </P>
<PRE>    A.sym_submatrix(f,l)            //   assumes fr=fc and lr=lc
    A.rows(f,l)                     //   select rows
    A.row(f)                        //   select single row
    A.columns(f,l)                  //   select columns
    A.column(f)                     //   select single column
</PRE>

<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>

<P>In each case <TT>f</TT> and <TT>l</TT> mean the first and last row or column
to be selected (starting at 1). </P>
<P>I allow <TT>l = f-1</TT> to indicate that a matrix of zero rows or columns
is to be returned. </P>
<P>If submatrix or its variant occurs on the right hand side of an <TT>=</TT>
or <TT>&lt;&lt;</TT> or within an expression think of its type as follows </P>
<PRE>    A.submatrix(fr,lr,fc,lc)           If A is RowVector or
                                       ColumnVector then same type
                                       otherwise type Matrix
    A.sym_submatrix(f,l)               Same type as A
    A.rows(f,l)                        Type Matrix
    A.row(f)                           Type RowVector
    A.columns(f,l)                     Type Matrix
    A.column(f)                        Type ColumnVector
</PRE>

<P>If submatrix or its variant appears on the left hand side of <TT>=</TT> or
<TT>&lt;&lt;</TT> , think of its type being Matrix. Thus <TT>L.row(1)</TT>
where <TT>L</TT> is LowerTriangularMatrix expects <TT>L.ncols()</TT> elements
even though it will use only one of them. If you are using <TT>=</TT> the
program will check for no loss of data.</P>

<P>A submatrix can appear on the left-hand side of <TT>+=</TT> or <TT>-=</TT>
with a matrix expression on the right-hand side. It can also appear on the
left-hand side of <TT>+=</TT>, <TT>-=</TT>, <TT>*=</TT> or <TT>/=</TT> with a
Real on the right-hand side. In each case there must be no loss of information.
The <a href="#binary">SP_eq function</a> is also defined for a submatrix.
</P>
<P>The <TT>row</TT> version can appear on the left hand side of
<TT>&lt;&lt;</TT> for <A HREF="#entering">loading literal data</A> into a row.
Load only the number of elements that are actually going to be stored in
memory. </P>
<P><b>Do not use the <TT>+=</TT> and <TT>-=</TT> operations with a submatrix of a
SymmetricMatrix or BandSymmetricMatrix on the LHS.</b></P>
<P>You can't pass a submatrix (or any of its variants) as a reference 
non-constant matrix in a function argument. For example, the following will not 
work:</P>
<pre>   void YourFunction(Matrix&amp; A);
   ...
   Matrix B(10,10);
   YourFunction(B.submatrix(1,5,1,5))    // won't compile</pre>
<P>If you are are using the submatrix facility to build a matrix from a small
number of components, consider instead using the <A
HREF="#binary">concatenation operators</A>. </P>
<H2><A NAME="dimen"></A>3.12 Change dimensions</H2>
<P CLASS="small"><A HREF="#ch_type">next</A> - <A HREF="#ch_type">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The following operations change the dimensions of a matrix. The values of the 
elements are lost, for <TT>resize</TT>. <tt>resize_keep</tt> keeps element 
values and zeros the elements in the matrix with the new size that are not in 
the matrix with the old size.</P>
<PRE>    A.resize(nrows,ncols);        // for type Matrix or nricMatrix
    A.resize(n);                  // for other types, except Band
    A.resize(n,lower,upper);      // for BandMatrix
    A.resize(n,lower);            // for LowerBandMatrix
    A.resize(n,upper);            // for UpperBandMatrix
    A.resize(n,lower);            // for SymmetricBandMatrix
    A.resize(B);                  // set dims to those of B
    A.cleanup();                  // resize to zero dimensions
    A.resize_keep(nrows,ncols);   // for type Matrix or nricMatrix, keep values
    A.resize_keep(n);             // for other types, except Band, keep values
</PRE>

<P>Use <TT>A.cleanup()</TT> to set the dimensions of <TT>A</TT> to zero and
release all the heap memory. </P>
<P><TT>A.resize(B)</TT> sets the dimensions of <TT>A</TT> to those of a matrix
<TT>B</TT>. This includes the band-width in the case of a band matrix. It is an
error for <TT>A</TT> to be a band matrix and <TT>B</TT> not a band matrix (or
diagonal matrix). </P>
<P>Remember that <TT>resize</TT> destroys values. If you want to resize, but 
keep the values in the bit that is left use <tt>resize_keep</tt>.</P>
<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>

<H2><A NAME="ch_type"></A>3.13 Change type</H2>
<P CLASS="small"><A HREF="#solve">next</A> - <A HREF="#solve">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The following functions interpret the elements of a matrix (stored row by
row) to be a vector or matrix of a different type. Actual copying is usually
avoided where these occur as part of a more complicated expression. </P>
<PRE>    A.as_row()
    A.as_column()
    A.as_diagonal()
    A.as_matrix(nrows,ncols)
    A.as_scalar()
</PRE>

<P>The expression <TT>A.as_scalar()</TT> is used to convert a 1 x 1 matrix to a
scalar.</P>
<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>
<H2><A NAME="solve"></A>3.14 Multiple matrix solve</H2>
<P CLASS="small"><A HREF="#memory">next</A> - <A HREF="#memory">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>To solve the matrix equation <TT>Ay = b</TT> where <TT>A</TT> is a square
matrix of equation coefficients, <TT>y</TT> is a column vector of values to be
solved for, and <TT>b</TT> is a column vector, use the code </P>
<PRE>    int n = something
    Matrix A(n,n); ColumnVector b(n);
    ... put values in A and b
    ColumnVector y = A.i() * b;       // solves matrix equation
</PRE>

<P>The following notes are for the case where you want to solve more than one
matrix equation with different values of <TT>b</TT> but the same <TT>A</TT>. Or
where you want to solve a matrix equation and also find the determinant of
<TT>A</TT>. In these cases you probably want to avoid repeating the LU
decomposition of <TT>A</TT> for each solve or determinant calculation. </P>
<P>If <TT>A</TT> is a square or symmetric matrix use </P>
<PRE>    ColumnVector p, q;    
    ...
    CroutMatrix X = A;                // carries out LU decomposition
    ColumnVector Ap = X.i()*p; ColumnVector Aq = X.i()*q;
    LogAndSign ld = X.log_determinant();
</PRE>

<P>rather than </P>
<PRE>    ColumnVector p, q;    
    ...
    ColumnVector Ap = A.i()*p; ColumnVector Aq = A.i()*q;
    LogAndSign ld = A.log_determinant();
</PRE>

<P>since each operation will repeat the LU decomposition. </P>
<P>If <TT>A</TT> is a BandMatrix or a SymmetricBandMatrix begin with </P>
<PRE>    BandLUMatrix X = A;               // carries out LU decomposition
</PRE>

<P>A CroutMatrix or BandLUMatrix can be copied and you can have a constructor 
with no parameters (use <tt>=</tt> to give it values). They work with
<a href="#memory"><i>release()</i>, <i>release_and_delete</i></a> and
<a href="#unspec">GenericMatrix</a> and <a href="#memory">ReturnMatrix</a>. You 
can't do any other manipulation apart from taking the inverse or solving with
<tt>i()</tt>, or finding the determinant or log determinant. See the <a href="#function">function summary list</a> 
for accessing the internals of a CroutMatrix or BandLUMatrix.</P>
<P>You can alternatively use </P>
<PRE>    LinearEquationSolver X = A;
</PRE>

<P>This will choose the most appropriate decomposition of <TT>A</TT>. That is,
the band form if <TT>A</TT> is banded; the Crout decomposition if <TT>A</TT> is
square or symmetric and no decomposition if <TT>A</TT> is triangular or 
diagonal. It doesn't know about positive definite matrices so won't use 
Cholesky.</P>
<H2><A NAME="memory"></A>3.15 Memory management</H2>
<P CLASS="small"><A HREF="#efficien">next</A> - <A HREF="#efficien">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The package does not support delayed copy. Several strategies are required
to prevent unnecessary matrix copies. </P>
<P>Where a matrix is called as a function argument use a constant reference.
For example </P>
<PRE>    YourFunction(const Matrix&amp; A)
</PRE>

<P>rather than </P>
<PRE>    YourFunction(Matrix A)
</PRE>

<P>Skip the rest of this section on your first reading. </P>
<P>A second place where it is desirable to avoid unnecessary copies is when a
function is returning a matrix. Matrices can be returned from a function with
the return command as you would expect. However these may incur one and
possibly two copyings of the matrix. To avoid this use the following
instructions. </P>
<P>Make your function of type ReturnMatrix . Then precede the return statement
with a <i>release</i> statement (or a <i>release_and_delete</i> statement if the matrix was
created with new). For example </P>
<PRE>    ReturnMatrix MakeAMatrix()
    {
       Matrix A;                // or any other matrix type
       ......
       A.release(); return A;
    }
</PRE>

<P>or </P>
<PRE>    ReturnMatrix MakeAMatrix()
    {
       Matrix* m = new Matrix;
       ......
       m-&gt;release_and_delete(); return *m;
    }
</PRE>

<P>If your compiler objects to this code, replace the return statements with 
</P>
<PRE>    return A.for_return();
</PRE>

<P>or </P>
<PRE>    return m-&gt;for_return();</PRE>

<HR>
<P align="center"><B>Do not forget to make the function of type ReturnMatrix.</B></P>
<HR>
<P>In particular, <b>don't</b> do</P>
<pre>    Matrix MakeAMatrix()
    {
       Matrix A;                // or any other matrix type
       ......
       A.release(); return A;   // will compile but could give <b>wrong answers</b>.
    }</pre>
<P>since with some compilers <tt>A</tt> might retain its <i>released</i> status 
after being returned.</P>
<P>You can also use <TT>.release()</TT> or <TT>-&gt;release_and_delete()</TT> to
allow a matrix expression to recycle space. Suppose you call </P>
<PRE>    A.release();
</PRE>

<P>just before <TT>A</TT> is used just once in an expression. Then the memory
used by <TT>A</TT> is either returned to the system or reused in the
expression. In either case, <TT>A</TT>'s memory is destroyed. This procedure
can be used to improve efficiency and reduce the use of memory. </P>
<P>Use <TT>-&gt;release_and_delete</TT> for matrices created by new if you want
to completely delete the matrix after it is accessed.</P>
<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>
<H2><A NAME="efficien"></A>3.16 Efficiency</H2>
<P CLASS="small"><A HREF="#output">next</A> - <A HREF="#output">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The package tends to be not very efficient for dealing with matrices with
short rows. This is because some administration is required for accessing rows
for a variety of types of matrices. To reduce the administration a special
multiply routine is used for rectangular matrices in place of the generic one.
Where operations can be done without reference to the individual rows (such as
adding matrices of the same type) appropriate routines are used. </P>
<P>When you are using small matrices (say smaller than 10 x 10) you may find it
faster to use rectangular matrices rather than the triangular or symmetric
ones. </P>
<H2><A NAME="output"></A>3.17 Output</H2>
<P CLASS="small"><A HREF="#unspec">next</A> - <A HREF="#unspec">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>To print a matrix use an expression like </P>
<PRE>   Matrix A, B C;
   ......
   cout &lt;&lt; setw(10) &lt;&lt; setprecision(5) &lt;&lt; A &lt;&lt; endl;
   cout &lt;&lt; setw(10) &lt;&lt; setprecision(5) &lt;&lt; scientific &lt;&lt; B &lt;&lt; endl;
   cout &lt;&lt; setw(10) &lt;&lt; setprecision(5) &lt;&lt; fixed &lt;&lt; C &lt;&lt; endl;
</PRE>

<P>This will work only with systems that support the standard input/output
routines including manipulators. The scientific and fixed manipulators won't 
work with Visual C++, version 6.</P>
<P>You need to #include the files iostream.h,
iomanip.h, newmatio.h in your C++ source files that use this facility. The
files iostream.h, iomanip.h will be included automatically if you include the
statement <TT>#define WANT_STREAM</TT> at the beginning of your source file. So
you can begin your file with either </P>
<PRE>   #define WANT_STREAM
   #include &quot;newmatio.h&quot;
</PRE>

<P>or </P>
<PRE>   #include &lt;iostream&gt;
   #include &lt;iomanip&gt;
   #include &quot;newmatio.h&quot;
</PRE>

<P>The present version of this routine is useful only for matrices small enough
to fit within a page or screen width. </P>
<P>To print several vectors or matrices in columns use a <A
HREF="#binary">concatenation operator</A>: </P>
<PRE>   ColumnVector A, B;
   .....
   cout &lt;&lt; setw(10) &lt;&lt; setprecision(5) &lt;&lt; (A | B) &lt;&lt; endl;
</PRE>

<H2><A NAME="unspec"></A>3.18 Unspecified type</H2>
<P CLASS="small"><A HREF="#cholesky">next</A> - <A HREF="#cholesky">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>Skip this section on your first reading. </P>
<P>If you want to work with a matrix of unknown type, say in a function. You
can construct a matrix of type <TT>GenericMatrix</TT>. Eg </P>
<PRE>   Matrix A;
   .....                                  // put some values in A
   GenericMatrix GM = A;
</PRE>

<P>A GenericMatrix matrix can be used anywhere where a matrix expression can be
used and also on the left hand side of an <TT>=</TT>. You can pass any type of 
matrix to a <TT>const
GenericMatrix&amp;</TT> argument in a function. However most scalar functions
including nrows(), ncols(), type() and element access do not work with it. Nor
does the ReturnMatrix construct. <tt><a href="#copy">swap</a></tt> does work with objects of type 
GenericMatrix. See also the paragraph on <A
HREF="#solve">LinearEquationSolver</A>. </P>
<P>An alternative and less flexible approach is to use BaseMatrix or
GeneralMatrix. </P>
<P>Suppose you wish to write a function which accesses a matrix of unknown type
including expressions (eg <TT>A*B</TT>). Then use a layout similar to the
following: </P>
<PRE>   void YourFunction(BaseMatrix&amp; X)
   {
      GeneralMatrix* gm = X.Evaluate();   // evaluate an expression
                                          // if necessary
      ........                            // operations on *gm
      gm-&gt;tDelete();                      // delete *gm if a temporary
   }
</PRE>

<P>See, as an example, the definitions of <TT>operator&lt;&lt;</TT> in
newmat9.cpp. </P>
<P>Under certain circumstances; particularly where <TT>X</TT> is to be used
just once in an expression you can leave out the <TT>Evaluate()</TT> statement
and the corresponding <TT>tDelete()</TT>. Just use <TT>X</TT> in the
expression. </P>
<P>If you know YourFunction will never have to handle a formula as its argument
you could also use </P>
<PRE>   void YourFunction(const GeneralMatrix&amp; X)
   {
      ........                            // operations on X
   }
</PRE>

<P>Do not try to construct a GeneralMatrix or BaseMatrix. </P>
<H2><A NAME="cholesky"></A>3.19 Cholesky
decomposition</H2>
<P CLASS="small"><A HREF="#qr">next</A> - <A HREF="#qr">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>Suppose <TT>S</TT> is symmetric and positive definite. Then there exists a
unique lower triangular matrix <TT>L</TT> such that <TT>L * L.t() = S</TT>. To
calculate this use </P>
<PRE>    SymmetricMatrix S;
    ......
    LowerTriangularMatrix L = Cholesky(S);
</PRE>

<P>If <TT>S</TT> is a symmetric band matrix then <TT>L</TT> is a band matrix
and an alternative procedure is provided for carrying out the decomposition: 
</P>
<PRE>    SymmetricBandMatrix S;
    ......
    LowerBandMatrix L = Cholesky(S);</PRE>

<p>See section <a href="#upd_chol">3.32</a> on updating a Cholesky 
decomposition.<br>
</p>

<H2><A NAME="qr"></A>3.20 QR decomposition</H2>
<P CLASS="small"><A HREF="#svd">next</A> - <A HREF="#svd">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>This is a variant on the usual QR transformation. </P>
<P>Start with matrix (dimensions shown to left and below the matrix)</P>
<PRE>       / 0    0 \      s
       \ X    Y /      n

         s    t
</PRE>

<P>Our version of the QR decomposition multiplies this matrix by an orthogonal
matrix Q to get </P>
<PRE>       / U    M \      s
       \ 0    Z /      n

         s    t
</PRE>

<P>where <TT>U</TT> is upper triangular (the R of the QR transform). That is 
</P>
<PRE>      Q  / 0   0 \  =  / U   M \
         \ X   Y /     \ 0   Z / </PRE>

<P>This is good for solving least squares problems: choose b (matrix or column
vector) to minimise the sum of the squares of the elements of </P>
<PRE>         Y - X*b
</PRE>

<P>Then choose <TT>b = U.i()*M;</TT> The residuals <TT>Y - X*b</TT> are in
<TT>Z</TT>. </P>
<P>This is the usual QR transformation applied to the matrix <TT>X</TT> with
the square zero matrix concatenated on top of it. It gives the same triangular
matrix as the QR transform applied directly to <TT>X</TT> and generally seems
to work in the same way as the usual QR transform. However it fits into the
matrix package better and also gives us the residuals directly. It turns out to
be essentially a modified Gram-Schmidt decomposition. </P>
<P>Two routines are provided in <I>newmat</I>: </P>
<PRE>    QRZ(X, U);
</PRE>

<P>replaces <TT>X</TT> by orthogonal columns and forms <TT>U</TT>. </P>
<PRE>    QRZ(X, Y, M);
</PRE>

<P>uses <TT>X</TT> from the first routine, replaces <TT>Y</TT> by <TT>Z</TT>
and forms <TT>M</TT>. There is also a routine</P>
<pre>    QRZ(X, Y, U, M);</pre>
<p>which does both of these.</p>
<P>To extend <TT>U</TT> to a square orthogonal matrix see the function for
<a href="#extend">extending an orthonormal set of columns</a>. </P>
<P>The are also routines <TT>QRZT(X, L)</TT>, <TT>QRZT(X, Y, M)</TT>
and <TT>QRZT(X, Y, L, M)</TT>
which do the same decomposition on the transposes of all these matrices. <tt>QRZT</tt>
replaces the routines <tt>HHDecompose</tt> in earlier versions of <i>newmat</i>. 
<tt>HHDecompose</tt> is
still defined but just calls <tt>QRZT</tt>.</P>
<P>For an example of the use of this decomposition see the file
<A HREF="#example">example.cpp</A>. </P>
<P>See the section on <a href="#upd_chol">updating a Cholesky decomposition</a> 
for updating <tt>U</tt>.</P>
<P>Alternatively to update <TT>U</TT> or <TT>L</TT> with a new block of data, <tt>X</tt>, one can use</P>
<pre>    updateQRZ(X, U);
</pre>
<p>or</p>
<pre>    updateQRZT(X, L);</pre>
<p>You can update the M matrix with</p>
<pre>    updateQRZ(X, Y, M);</pre>
<p>At present, there is no corresponding function for <tt>updateQRZT</tt>. Also 
note, at present, <tt>updateQRZ(X, Y, M)</tt> is not optimised for accessing
<a href="#mem_man2">contiguous locations</a>.</p>
<p>If you have carried out QRZ transforms on two blocks of data, <tt>X1</tt>, <tt>X2</tt>, 
(same number of columns in each) with the data to be fitted in <tt>Y1</tt>, <tt>
Y2</tt>, you can combine them as follows</p>
<pre>    // QRZ transforms of two blocks of data
    UpperTriangularMatrix U1, U2; Matrix M1, M2;
    QRZ(X1, U1); QRZ(X2, U2); QRZ(X1, Y1, M1); QRZ(X2, Y2, M2);

    // Now combine the results
    updateQRZ(U1, U2); updateQRZ(U1, M1, M2);</pre>
<p>Results are in <tt>U2</tt> and <tt>M2</tt>. The values in <tt>U1</tt> and <tt>
M1</tt> are modified so use copies if you still need the original values.</p>
<p>The functions <tt>updateQRZ(U1, U2)</tt> and <tt>updateQRZ(U1, M1, M2)</tt> 
are not optimised for accessing <a href="#mem_man2">contiguous locations</a>. 
This is probably not an issue, since usually all the data will fit into cache 
memory.&nbsp; </p>
<p><b>Notes on updateQRZ and updateQRZT:</b></p>
<ul>
  <li>Use these routines if your data is arriving in blocks or there is too much to 
fit into memory.</li>
  <li>Use <TT>QRZ</TT> or <TT>QRZT</TT> on the first block of data or you can use <TT>
  updateQRZ</TT> or 
<TT>updateQRZT</TT> if <TT>U</TT> or <TT>L</TT> is correctly dimensioned and set to zero.</li>
  <li>The block of data, <tt>X</tt>, will be modified by these routines but the 
modified data is probably not useful for further analysis.</li>
  <li>The signs of some of the rows of <TT>U</TT> or columns of <TT>L</TT> maybe the reverse of what you would 
expect (i.e. you get the negative of what you get from QRZ or QRZT applied to 
	the whole dataset). Usually this will not be a problem.</li>
	<li>The second version of updateQRZ is sometimes a more useful way of 
	combining blocks, especially if you want to use the <i>reduce</i> function 
	in <i>Intel's Threaded Building Blocks</i> to parallelise QRZ. (Contact me 
	for what you need to do to make Newmat compatible with <i>Threaded Building 
	Blocks</i>).</li>
  <li>See the <a href="#function">function summary list</a> for the older 
depreciated function names.</li>
</ul>
<p>&nbsp;</p>
<H2><A NAME="svd"></A>3.21 Singular value
decomposition</H2>
<P CLASS="small"><A HREF="#evalues">next</A> - <A HREF="#evalues">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The singular value decomposition of an <I>m</I> x <I>n</I> <TT>Matrix</TT>
<TT>A</TT> (where <I>m</I> &gt;= <I>n</I>) is a decomposition </P>
<PRE>    A  = U * D * V.t()
</PRE>

<P>where <TT>U</TT> is <I>m</I> x <I>n</I> with <TT>U.t() * U</TT> equalling
the identity, <TT>D</TT> is an <I>n</I> x <I>n </I><TT>DiagonalMatrix</TT> and
<TT>V</TT> is an <I>n</I> x <I>n</I> orthogonal matrix (type <TT>Matrix</TT> in
<I>Newmat</I>). </P>
<P>Singular value decompositions are useful for understanding the structure of
ill-conditioned matrices, solving least squares problems, and for finding the
eigenvalues of <TT>A.t() * A</TT>. </P>
<P>To calculate the singular value decomposition of <TT>A</TT> (with <I>m</I>
&gt;= <I>n</I>) use one of </P>
<PRE>    SVD(A, D, U, V);                  // U = A is OK
    SVD(A, D);
    SVD(A, D, U);                     // U = A is OK
    SVD(A, D, U, false);              // U (can = A) for workspace only
    SVD(A, D, U, V, false);           // U (can = A) for workspace only
</PRE>

<P>where <TT>A</TT>, <TT>U</TT> and <TT>V</TT> are of type <TT>Matrix</TT> and
<TT>D</TT> is a <TT>DiagonalMatrix</TT>. The values of <TT>A</TT> are not
changed unless <TT>A</TT> is also inserted as the third argument.</P>

<P>The elements of <tt>D</tt> are sorted in <i>descending</i> order.</P>

<P>To extend <TT>U</TT> to a square orthogonal matrix see the function for
<a href="#extend">extending an orthonormal set of columns</a>. </P>
<P CLASS="small">Remember that the SVD decomposition is not completely unique. The signs of the elements in a column of <TT>U</TT> may be reversed
if the signs in the corresponding column in <TT>V</TT> are reversed. If a
number of the singular values are identical one can apply an orthogonal
transformation to the corresponding columns of <TT>U</TT> and the corresponding
columns of <TT>V</TT>.</P>
<P CLASS="small">If <i>m</i> &lt; <i>n</i> apply the SVD transform to the transpose 
of A and swap U and V. If necessary, extend U to a square matrix as described 
above.</P>
<H2><A NAME="evalues"></A>3.22 Eigenvalue
decomposition</H2>
<P CLASS="small"><A HREF="#sorting">next</A> - <A HREF="#sorting">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>An eigenvalue decomposition of a SymmetricMatrix <TT>A</TT> is a
decomposition </P>
<PRE>    A  = V * D * V.t()
</PRE>

<P>where <TT>V</TT> is an orthogonal matrix (type <TT>Matrix</TT> in
<I>Newmat</I>) and <TT>D</TT> is a DiagonalMatrix. </P>
<P>Eigenvalue analyses are used in a wide variety of engineering, statistical
and other mathematical analyses. </P>
<P>The package includes two algorithms: Jacobi and Householder. The first is
extremely reliable but much slower than the second. </P>
<P>The code is adapted from routines in <I>Handbook for Automatic Computation,
Vol II, Linear Algebra</I> by Wilkinson and Reinsch, published by Springer
Verlag. </P>
<PRE>    Jacobi(A,D,S,V);                  // A, S symmetric; S is workspace,
                                      //    S = A is OK; V is a matrix
    Jacobi(A,D);                      // A symmetric
    Jacobi(A,D,S);                    // A, S symmetric; S is workspace,
                                      //    S = A is OK
    Jacobi(A,D,V);                    // A symmetric; V is a matrix

    eigenvalues(A,D);                 // A symmetric
    eigenvalues(A,D,S);               // A, S symmetric; S is for back
                                      //    transforming, S = A is OK
    eigenvalues(A,D,V);               // A symmetric; V is a matrix
</PRE>

<P>where <TT>A</TT>, <TT>S</TT> are of type <TT>SymmetricMatrix</TT>,
<TT>D</TT> is of type <TT>DiagonalMatrix</TT> and <TT>V</TT> is of type
<TT>Matrix</TT>. The values of <TT>A</TT> are not changed unless <TT>A</TT> is
also inserted as the third argument. If you need eigenvectors use one of the
forms with matrix <TT>V</TT>. The eigenvectors are returned as the columns of
<TT>V</TT>.</P>

<P>The elements of <tt>D</tt> are sorted in <i>ascending</i> order.</P>
<P CLASS="small">Remember that an eigenvalue decomposition is not completely
unique - see the comments about the <A HREF="#svd">SVD</A>
decomposition.</P>
<P CLASS="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names. </P>
<H2><A NAME="sorting"></A>3.23 Sorting</H2>
<P CLASS="small"><A HREF="#fft">next</A> - <A HREF="#fft">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>To sort the values in a matrix or vector, <TT>A</TT>, (in general this
operation makes sense only for vectors and diagonal matrices) use one of</P>
<PRE>    sort_ascending(A);

    sort_descending(A);
</PRE>

<P>I use the quicksort algorithm. The algorithm is similar to that in
Sedgewick's algorithms in C++. If the sort seems to be failing (as quicksort
can do) an exception is thrown. </P>
<P>You will get incorrect results if you try to sort a band matrix - but why
would you want to sort a band matrix? </P>
<P class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</P>
<H2><A NAME="fft"></A>3.24 Fast Fourier transform 
</H2>
<P CLASS="small"><A HREF="#trigtran">next</A> - <A HREF="#trigtran">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<PRE>   FFT(X, Y, F, G);                         // F=X and G=Y are OK
</PRE>

<P>where <TT>X</TT>, <TT>Y</TT>, <TT>F</TT>, <TT>G</TT> are ColumnVectors.
<TT>X</TT> and <TT>Y</TT> are the real and imaginary input vectors; <TT>F</TT>
and <TT>G</TT> are the real and imaginary output vectors. The lengths of
<TT>X</TT> and <TT>Y</TT> must be equal and should be the product of numbers
less than about 10 for fast execution. </P>
<P>The formula is </P>
<PRE>          n-1
   h[k] = SUM  z[j] exp (-2 pi i jk/n)
          j=0
</PRE>

<P>where <TT>z[j]</TT> is stored complex and stored in <TT>X(j+1)</TT> and
<TT>Y(j+1)</TT>. Likewise <TT>h[k]</TT> is complex and stored in
<TT>F(k+1)</TT> and <TT>G(k+1)</TT>. The fast Fourier algorithm takes order <I>
n</I> log(<I>n</I>) operations (for <I>good</I> values of <I>n</I>) rather than
<I>n</I>**2 that straight evaluation (see the file <TT>tmtf.cpp</TT>) takes. 
</P>
<P>I use one of two methods: </P>
<UL>
<LI>A program originally written by Sande and Gentleman. This requires that
<I>n</I> can be expressed as a product of small numbers.</LI>
<LI>A method of Carl de Boor (1980), <I>Siam J Sci Stat Comput</I>, pp 173-8.
The sines and cosines are calculated explicitly. This gives better accuracy, at
an expense of being a little slower than is otherwise possible. This is slower
than the Sande-Gentleman program but will work for all <I>n</I> --- although it
will be very slow for <I>bad</I> values of <I>n</I>.</LI>
</UL>
<P>Related functions </P>
<PRE>   FFTI(F, G, X, Y);                        // X=F and Y=G are OK
   RealFFT(X, F, G);
   RealFFTI(F, G, X);
</PRE>

<P><TT>FFTI</TT> is the inverse transform for <TT>FFT</TT>. <TT>RealFFT</TT> is
for the case when the input vector is real, that is <TT>Y = 0</TT>. I assume
the length of <TT>X</TT>, denoted by <I>n</I>, is <I>even</I>. That is,
<I>n</I> must be divisible by 2. The program sets the lengths of <TT>F</TT> and
<TT>G</TT> to <I>n</I>/2 + 1. <TT>RealFFTI</TT> is the inverse of
<TT>RealFFT</TT>. </P>
<P>See also the section on fast <A HREF="#trigtran">trigonometric
transforms</A>.</P>
<P>There  are also two dimensional versions</P>
<pre>   FFT2(X, Y, F, G);                       // F=X and G=Y are OK
   FFT2I(F, G, X, Y);                      // inverse, X=F and Y=G are OK</pre>
<p>where <TT>X</TT>, <TT>Y</TT>, <TT>F</TT>, <TT>G</TT> are of type Matrix.
<TT>X</TT> and <TT>Y</TT> are the real and imaginary input matrices; <TT>F</TT>
and <TT>G</TT> are the real and imaginary output matrices. The dimensions of <TT>Y</TT> 
must be the same as those of <TT>X</TT> and should be the product of numbers
less than about 10 for fast execution. </p>
<P>The formula is </P>
<PRE>            m-1 n-1
   h[p,q] = SUM SUM z[j,k] exp (-2 pi i (jp/m + kq/n))
            j=0 k=0
</PRE>

<P>where <TT>z[j,k]</TT> is stored complex and stored in <TT>X(j+1,k+1)</TT> and
<TT>Y(j+1,k+1)</TT> and <TT>X</TT> and <TT>Y</TT> have dimension <TT>m x n</TT>. Likewise <TT>h[p,q]</TT> is complex and stored in
<TT>F(p+1,q+1)</TT> and <TT>G(p+1,q+1)</TT>. 
</P>
<P>&nbsp;</P>
<H2><A NAME="trigtran"></A>3.25 Fast trigonometric
transforms</H2>
<P CLASS="small"><A HREF="#nric">next</A> - <A HREF="#nric">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>These are the sin and cosine transforms as defined by Charles Van Loan
(1992) in <I>Computational frameworks for the fast Fourier transform</I>
published by SIAM. See page 229. Some other authors use slightly different
conventions. All the functions call the <A HREF="#fft">fast Fourier
transforms</A> and require an <I> even</I> transform length, denoted by
<I>m</I> in these notes. That is, <I>m</I> must be divisible by 2. As with the
FFT <I>m</I> should be the product of numbers less than about 10 for fast
execution. </P>
<P>The functions I define are </P>
<PRE>   DCT(U,V);                // U, V are ColumnVectors, length <I>m+1</I>
   DCT_inverse(V,U);        // inverse of DCT
   DST(U,V);                // U, V are ColumnVectors, length <I>m+1</I>
   DST_inverse(V,U);        // inverse of DST
   DCT_II(U,V);             // U, V are ColumnVectors, length <I>m</I>
   DCT_II_inverse(V,U);     // inverse of DCT_II
   DST_II(U,V);             // U, V are ColumnVectors, length <I>m</I>
   DST_II_inverse(V,U);     // inverse of DST_II
</PRE>

<P>where the first argument is the input and the second argument is the output.
<TT>V = U</TT> is OK. The length of the output ColumnVector is set by the
functions. </P>
<P>Here are the formulae: </P>
<H3>DCT</H3>
<PRE>                   m-1                             k
   v[k] = u[0]/2 + SUM { u[j] cos (pi jk/m) } + (-) u[m]/2
                   j=1
</PRE>

<P>for <TT>k = 0...m</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored in
<TT>U(j+1)</TT> and <TT>V(k+1)</TT>. </P>
<H3>DST</H3>
<PRE>          m-1
   v[k] = SUM { u[j] sin (pi jk/m) }
          j=1
</PRE>

<P>for <TT>k = 1...(m-1)</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored
in <TT>U(j+1)</TT> and <TT>V(k+1)</TT>and where <TT>u[0]</TT> and <TT>u[m]</TT>
are ignored and <TT>v[0]</TT> and <TT>v[m]</TT> are set to zero. For the
inverse function <TT>v[0]</TT> and <TT>v[m]</TT> are ignored and <TT>u[0]</TT>
and <TT>u[m]</TT> are set to zero. </P>
<H3>DCT_II</H3>
<PRE>          m-1
   v[k] = SUM { u[j] cos (pi (j+1/2)k/m) }
          j=0
</PRE>

<P>for <TT>k = 0...(m-1)</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored
in <TT>U(j+1)</TT> and <TT>V(k+1)</TT>. </P>
<H3>DST_II</H3>
<PRE>           m
   v[k] = SUM { u[j] sin (pi (j-1/2)k/m) }
          j=1
</PRE>

<P>for <TT>k = 1...m</TT>, where <TT>u[j]</TT> and <TT>v[k]</TT> are stored in
<TT>U(j)</TT> and <TT>V(k)</TT>. </P>
<P>Note that the relationship between the subscripts in the formulae and those
used in <I>newmat</I> is different for DST_II (and DST_II_inverse). </P>
<H2><A NAME="nric"></A>3.26 Interface to Numerical
Recipes in C</H2>
<P CLASS="small"><A HREF="#except">next</A> - <A HREF="#except">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>This section refers to <i>Numerical Recipes in C</i>. This section is <b>not</b> 
relevant to <i>Numerical Recipes in C++</i>. I'll put a note on the website soon 
on how to interface with <i>Numerical Recipes in C++.</i></P>
<P>This package can be used with the vectors and matrices defined in
<I>Numerical Recipes in C</I>. You need to edit the routines in Numerical
Recipes so that the elements are of the same type as used in this package. Eg
replace float by double, vector by dvector and matrix by dmatrix, etc. You may
need to edit the function definitions to use the version acceptable to your
compiler (if you are using the first edition of NRIC). You may need to enclose
the code from Numerical Recipes in <TT>extern &quot;C&quot; { ... }</TT>. You
will also need to include the matrix and vector utility routines. </P>
<P>Then any vector in Numerical Recipes with subscripts starting from 1 in a
function call can be accessed by a RowVector, ColumnVector or DiagonalMatrix in
the present package. Similarly any matrix with subscripts starting from 1 can
be accessed by an nricMatrix in the present package. The class nricMatrix is
derived from Matrix and can be used in place of Matrix. In each case, if you
wish to refer to a RowVector, ColumnVector, DiagonalMatrix or nricMatrix
<TT>X</TT> in an function from Numerical Recipes, use <TT>X.nric()</TT> in the
function call. </P>
<P>Numerical Recipes cannot change the dimensions of a matrix or vector. So
matrices or vectors must be correctly dimensioned before a Numerical Recipes
routine is called. </P>
<P>For example </P>
<PRE>   SymmetricMatrix B(44);
   .....                             // load values into B
   nricMatrix BX = B;                // copy values to an nricMatrix
   DiagonalMatrix D(44);             // Matrices for output
   nricMatrix V(44,44);              //    correctly dimensioned
   int nrot;
   jacobi(BX.nric(),44,D.nric(),V.nric(),&amp;nrot);
                                     // jacobi from NRIC
   cout &lt;&lt; D;                        // print eigenvalues
</PRE>

<H2><A NAME="except"></A>3.27 Exceptions</H2>
<P CLASS="small"><A HREF="#cleanup">next</A> - <A HREF="#cleanup">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>Here is the class structure for exceptions: </P>
<PRE>Exception
  Logic_error
    ProgramException                 miscellaneous matrix error
    IndexException                   index out of bounds
    VectorException                  unable to convert matrix to vector
    NotSquareException               matrix is not square (invert, solve)
    SubMatrixDimensionException      out of bounds index of submatrix
    IncompatibleDimensionsException  (multiply, add etc)
    NotDefinedException              operation not defined (eg &lt;)
    CannotBuildException             copying a matrix where copy is undefined
    InternalException                probably an error in newmat
  Runtime_error
    NPDException                     matrix not positive definite (Cholesky)
    ConvergenceException             no convergence (e-values, non-linear, sort)
    SingularException                matrix is singular (invert, solve)
    SolutionException                no convergence in solution routine
    OverflowException                floating point overflow
  Bad_alloc                          out of space (new fails)
</PRE>

<P>I have attempted to mimic the exception class structure in the C++ standard
library, by defining the Logic_error and Runtime_error classes.</P>
<P>Suppose you have edited <tt>include.h</tt> to use my <i>simulated</i>
exceptions or to <i>disable</i> exceptions. If there is no catch statement or exceptions are disabled then my
<TT>Terminate()</TT> function in <TT>myexcept.h</TT> is called when you throw an
exception. This prints out
an error message, the dimensions and types of the matrices involved, the name
of the routine detecting the exception, and any other information set by the
<A HREF="#error">Tracer</A> class. Also see the section on <A
HREF="#error">error messages</A> for additional notes on the messages generated
by the exceptions. </P>
<P>You can also print this information in a <i>catch</i> clause by printing <TT>Exception::what()</TT>. 
</P>
<P>If you are using <i> compiler supported</i> exceptions then see the section
on <a href="#except_1">catching exceptions</a>.&nbsp; 
</P>
<P>See the file <TT>test_exc.cpp</TT> as an example of catching an exception
and printing the error message. </P>
<P>The 08 version of newmat defined a member function <TT>void
SetAction(int)</TT> to help customise the action when an exception is called.
This has been deleted in the 09 and later versions. Now include an instruction
such as <TT>cout &lt;&lt; Exception::what() &lt;&lt; endl;</TT> in the
<TT>Catch</TT> or <TT>CatchAll</TT> block to determine the action. </P>
<P>The library includes the alternatives of using the inbuilt exceptions
provided by a compiler, simulating exceptions, or disabling exceptions. See
<A HREF="#custom">customising</A> for selecting the correct exception option. 
</P>
<P>The rest of this section describes my partial simulation of exceptions for
compilers which do not support C++ exceptions. Skip the rest of this section and 
the next section if you are using compiler supported exceptions. I use Carlos Vidal's article in
the September 1992 <I>C Users Journal</I> as a starting point. </P>
<P>Newmat does a partial clean up of memory following throwing an exception -
see the next section. However, the present version will leave a little heap
memory unrecovered under some circumstances. I would not expect this to be a
major problem, but it is something that needs to be sorted out. </P>
<P>The functions/macros I define are Try, Throw, Catch, CatchAll and
CatchAndThrow. Try, Throw, Catch and CatchAll correspond to try, throw, catch
and catch(...) in the C++ standard. A list of Catch clauses must be terminated
by either CatchAll or CatchAndThrow but not both. Throw takes an Exception as
an argument or takes no argument (for passing on an exception). I do not have a
version of Throw for specifying which exceptions a function might throw. Catch
takes an exception class name as an argument; CatchAll and CatchAndThrow don't
have any arguments. Try, Catch and CatchAll must be followed by blocks enclosed
in curly brackets. </P>
<P>I have added another macro ReThrow to mean a rethrow, Throw(). This was
necessary to enable the package to be compatible with both my exception package
and C++ exceptions. </P>
<P>If you want to throw an exception, use a statement like </P>
<PRE>   Throw(Exception(&quot;Error message\n&quot;));
</PRE>

<P>It is important to have the exception declaration in the Throw statement,
rather than as a separate statement. </P>
<P>All exception classes must be derived from the class, Exception, defined in
newmat and can contain only static variables. See the examples in newmat if you
want to define additional exceptions.</P>
<P>Note that the simulation exception mechanism does not work if you define
arrays of matrices. </P>
<H2><A NAME="cleanup"></A>3.28 Cleanup after an
exception</H2>
<P CLASS="small"><A HREF="#nonlin">next</A> - <A HREF="#nonlin">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>This section is about the <i>simulated exceptions</i> used in newmat. It is
<b>irrelevant</b> if you are using the exceptions built into a compiler or have set
the disable-exceptions option. </P>
<P>The simulated exception mechanisms in newmat are based on the C functions
setjmp and longjmp. These functions do not call destructors so can lead to
garbage being left on the heap. (I refer to memory allocated by <I>new</I> as
heap memory). For example, when you call </P>
<PRE>   Matrix A(20,30);
</PRE>

<P>a small amount of space is used on the stack containing the row and column
dimensions of the matrix and 600 doubles are allocated on the heap for the
actual values of the matrix. At the end of the block in which A is declared,
the destructor for A is called and the 600 doubles are freed. The locations on
the stack are freed as part of the normal operations of the stack. If you leave
the block using a longjmp command those 600 doubles will not be freed and will
occupy space until the program terminates. </P>
<P>To overcome this problem newmat keeps a list of all the currently declared
matrices and its exception mechanism will return heap memory when you do a
Throw and Catch. </P>
<P>However it will not return heap memory from objects from other packages. 
</P>
<P>If you want the mechanism to work with another class you will have to do
four things: </P>
<OL>
<LI>derive your class from class Janitor defined in except.h; </LI>
<LI>define a function <TT>void CleanUp()</TT> in that class to return all heap
memory; </LI>
<LI>include the following lines in the class definition <PRE>      public:
         void* operator new(size_t size)
         { do_not_link=true; void* t = ::operator new(size); return t; }
         void operator delete(void* t) { ::operator delete(t); }
</PRE>

</LI>
<LI>be sure to include a copy constructor in you class definition, that is,
something like <PRE>      X(const X&amp;);
</PRE>

</LI>
</OL>
<P>Use <TT>CleanUp()</TT> rather than <TT>cleanup()</TT> since this is what is 
defined in class Janitor.</P>
<P>Note that the function <TT>CleanUp()</TT> does somewhat the same duties as
the destructor. However <TT>CleanUp()</TT> has to do the <I>cleaning</I> for
the class you are working with and also the classes it is derived from. So it
will often be wrong to use exactly the same code for both <TT>CleanUp()</TT>
and the destructor or to define your destructor as a call to
<TT>CleanUp()</TT>. </P>
<H2><A NAME="nonlin"></A>3.29 Non-linear
applications</H2>
<P CLASS="small"><A HREF="#stl">next</A> - <A HREF="#stl">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>Files solution.h, solution.cpp contain a class for solving for <I>x</I> in
<I>y</I> = <I>f</I>(<I>x)</I> where <I>x</I> is a one-dimensional continuous
monotonic function. This is not a matrix thing at all but is included because
it is a useful program and because it is a simpler version of the coding technique used
in the non-linear least squares. </P>
<P>Files newmatnl.h, newmatnl.cpp contain a series of classes for non-linear
least squares and maximum likelihood. These classes work on very well-behaved
functions but need upgrading for less well-behaved functions. I haven't followed 
the usual practice of inflating the values of the diagonal elements of the 
matrix of second derivatives. I originally thought I could avoid this if my 
program had a good line search. But this was wrong and when I use this program 
on all but the most well-behaved problems I run the fit first with the diagonal 
elements inflated by a factor of 2 to 5 and the critical value for the stopping 
criterion set to something like 50. Then rerun with with no inflation factor and 
critical value 0.0001.</P>
<P>Documentation for both of these is in the definition files. Simple examples
are in sl_ex.cpp, nl_ex.cpp and garch.cpp.</P>
<H2><A NAME="stl"></A>3.30 Standard template
library</H2>
<P CLASS="small"><A HREF="#namesp">next</A> - <A HREF="#namesp">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P>The standard template library (STL) is the set of <I>container templates</I>
(vector, deque, list etc) defined by the C++ standards committee. Newmat is
intended to be compatible with the STL in the sense that you can store matrices
in the standard containers. I have defined <A HREF="#binary"><TT>==</TT> and
inequality </A>operators which seem to be required by some versions of the STL.</P>
<P>If you want to use the container classes with Newmat please note </P>
<UL>
<LI>Don't use simulated exceptions. </LI>
<LI>Make sure the option <A HREF="#custom">DO_FREE_CHECK</A> is <EM>not</EM>
turned on. </LI>
<LI>You can store only one type of matrix in a container. If you want to use a
variety of types use the GenericMatrix type or store pointers to the matrices. 
</LI>
<LI>The vector and deque container templates like to copy their elements. For
the vector container this happens when you insert an element anywhere except at
the end or when you append an element and the current vector storage overflows.
Since Newmat does not have <I>copy-on-write</I> this could get very inefficient. </LI>
<LI>You won't be able to sort the container or do anything that would call an
inequality operator. </LI>
</UL>
<H2><A NAME="namesp"></A>3.31 Namespace</H2>
<P CLASS="small"><A HREF="#upd_chol">next</A> - <A HREF="#upd_chol">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></P>
<P><I>Namespace</I> is used to avoid name
clashes between different libraries. I have included the namespace capability.
Activate the line <TT>#define use_namespace</TT> in <TT>include.h</TT>. Then
include either the statement </P>
<PRE>   using namespace NEWMAT;
</PRE>

<P>at the beginning of any file that needs to access the newmat library or </P>
<PRE>   using namespace RBD_LIBRARIES;
</PRE>

<P>at the beginning of any file that needs to access all my libraries. </P>
<P>This works correctly with <A HREF="#borland">Borland</A> C++ version 5 and 
Builder 5 and 6. </P>
<P><A HREF="#microso">Microsoft</A> Visual C++ version 5 works in my example
and test files, but fails with apparently insignificant changes (it may be more
reliable if you have applied service pack 3). If you #include
&quot;newmatap.h&quot;, but no other newmat include file, then also #include
&quot;newmatio.h&quot;. It seems to work with <A HREF="#microso">Microsoft</A>
Visual C++ version 6 <I>if</I> you have applied at least service pack 2.</P>
<P>My use of namespace does not work with <A HREF="#gcc">Gnu g++</A> version 
2.8.1 but does work with versions 3.x. </P>
<P>I have defined the following namespaces: </P>
<UL>
<LI>RBD_COMMON for functions and classes used by most of my libraries </LI>
<LI>NEWMAT for the newmat library </LI>
<LI>RBD_LIBRARIES for all my libraries </LI>
</UL>
<H2><A NAME="upd_chol"></A>3.32 Updating the Cholesky matrix</H2>


<p CLASS="small"><A HREF="#RealStarStar">next</A> - <A HREF="#RealStarStar">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></p>


<p>Suppose <tt>X</tt> is matrix and <tt>U</tt> has been formed with either</p>


<pre>   SymmetricMatrix A; A &lt;&lt; X.t() * X;
   UpperTriangularMatrix U = Cholesky(A).t();</pre>
<p>or</p>
<pre>   UpperTriangularMatrix U;
   QRZ(X, U);</pre>
<p>See sections <a href="#cholesky">3.19</a> and <a href="#qr">3.20</a>.</p>
<p>Now suppose we want to append an extra row to <tt>X</tt> or delete a row from <tt>X</tt> 
or rearrange the columns of <tt>X.</tt> The functions described here allow you to 
update <tt>U</tt> without&nbsp;&nbsp;recalculating it.</p>


<pre>   update_Cholesky(U, x);                    // x is a RowVector
   downdate_Cholesky(U, x);                  // x is a RowVector
   right_circular_update_Cholesky(U, j, k);  // j and k are integers
   left_circular_update_Cholesky(U, j, k);   // j and k are integers</pre>
<p><tt>update_Cholesky</tt> carries out the modification of <tt>U</tt> corresponding to appending 
an extra row <tt>x</tt> to <tt>X</tt>.</p>


<p> <tt>downdate_Cholesky</tt> carries out the modification corresponding 
to removing a row <tt>x</tt> from <tt>X</tt>. A <tt>ProgramException</tt> exception is 
thrown if the modification fails.</p>


<p><tt>right_circular_update_Cholesky</tt> supposes that columns <tt>j,j+1,...k</tt> of <tt>X</tt> are 
replaced by columns <tt>k,j,j+1,...k-1</tt>.</p>


<p><tt>left_circular_update_Cholesky</tt> supposes that columns <tt>j,j+1,...k</tt> of <tt>X</tt> are replaced 
by columns <tt>j+1,...k,j</tt>.</p>


<p>These functions are based on a contribution from Nick Bennett of 
Schlumberger-Doll Research. See also the LINPACK User's Guide, Chapter 10, 
Dongarra et. al., SIAM, Philadelphia, 1979.</p>


<p>Where you want to append a number of new rows consider using the update 
routine in section <a href="#qr">3.20</a>.</p>


<p class="small">See the <a href="#function">function summary list</a> for the older 
depreciated function names.</p>


<h2><a name="RealStarStar"></a>3.33 Accessing C functions</h2>


<p CLASS="small"><A HREF="#SimpleIntArray">next</A> - <A HREF="#SimpleIntArray">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></p>


<p>You have a <i>C</i> function that uses one and two dimensional arrays as 
vectors and matrices. You want to access it from <i>Newmat</i>.</p>


<p>One dimensional arrays are easy. Set up a <i>ColumnVector</i>, <i>RowVector</i> 
or <i>DiagonalMatrix</i> with the correct dimension and where the function has a
<tt>double*</tt> argument enter <tt>X.data()</tt> where <tt>X</tt> denotes the <i>
ColumnVector</i>, <i>RowVector</i> or <i>DiagonalMatrix</i>. (I am assuming you 
have left <i>Real</i> being <i>typedef</i>ed as a <i>double</i>). If you have a
<tt>const double*</tt> argument use <tt>X.const_data()</tt>.</p>


<p>You can't do this with two dimensional arrays where you have a <tt>double**</tt> 
argument. <i>Newmat</i> includes  classes <i>RealStarStar</i> and <i>ConstRealStarStar</i> for this situation. 
To access a <tt>Matrix A</tt> from a function <tt>c_function(double** a)</tt> 
use either</p>


<pre>   c_function(RealStarStar(A));</pre>
<p>or</p>
<pre>   RealStarStar a(A);
   c_function(a);</pre>
<p>If the argument is <tt>const double**</tt> use <i>ConstRealStarStar</i>.</p>


<p>The <tt>Matrix A</tt> must be&nbsp;correctly dimensioned and must not be 
resized or set equal to another matrix between setting up the <i>RealStarStar</i> 
object and its use in the function. Also the following construction will not 
work</p>


<pre>   double** a = RealStarStar(A);      // wrong
   c_function(a);</pre>
<p>since the <i>RealStarStar</i> structure will have been destroyed before you 
get to the second line.</p>


<h2><a name="SimpleIntArray"></a>3.34 Simple integer array class</h2>


<p CLASS="small"><a href="#extend">next</a> - <a href="#extend">skip</a> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></p>


<p>This is primarily for use within <i>Newmat</i>. You can set up a simple array 
of integers with the <i>SimpleIntArray</i> class. Here are the descriptions of 
the constructors and functions.</p>


<TABLE WIDTH=100%>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>SimpleIntArray A;</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Constructs int array of length zero</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>SimpleIntArray A(n);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Constructs int array of length <i>n</i> - 
individual elements are <i>not</i> initialised</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A = i;</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> sets values to <i>i</i> where <i>i</i> is an int 
variable</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A = B;</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> sets values to those of <i>B</i> where 
<i>B</i> is 
a SimpleIntArray, change size if necessary.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>n = A.size();</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> return the length of <i>A</i></TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A.resize(n);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> change the length of <i>A</i>, don't keep 
values</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A.resize_keep(n);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> change the length of <i>A</i>, do keep 
values; if length is being increased set new elements to zero.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>int x = A[i];</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> element access; <i>i</i> runs from 0 to 
<i>n</i>-1</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>int* d = A.data();</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> get beginning of data array</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>const int* d = A.data();</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> get beginning of data array as const 
int* if <i>A</i> is const</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>const int* d = A.const_data();</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> get beginning of data array as const 
int*</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A.cleanup()</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> resize to zero length</TD>
</TR>
</TABLE>

&nbsp;<h2><a name="extend"></a>3.35

Extend orthonormal set of columns</h2>
<p class = "small"><A HREF="#misc_fn">next</A> - <A HREF="#misc_fn">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></p>
<p>Suppose a <i>Matrix</i> <tt>A</tt>'s first <tt>n</tt> columns are orthonormal so that 
<tt>A.Columns(1,n).t() * A.Columns(1,n)</tt> is the identity matrix. Suppose
we want to fill out the remaining columns of <tt>A</tt> to make them orthonormal so 
that <tt>A.t() * A</tt> is the identity matrix. Then use the function</p>
<pre>   extend_orthonormal(A, n);</pre>
<p>Matrix <tt>A</tt> is then <i>replaced</i> by the matrix with the additional 
columns.</p>
<p>Use this function to extend <tt>U</tt> from the <a href="#qr">QRZ</a> or
<a href="#svd">SVD</a> decompositions to form a square (orthogonal) matrix.</p>
<p>Notes:</p>
<ul>
  <li>the first <tt>n</tt> columns of <tt>A</tt> must be orthonormal</li>
  <li><tt>n</tt> must be less or equal to the number of columns of <tt>A</tt></li>
  <li>the number of columns of <tt>A</tt> must be less than or equal to the number of 
  rows of <tt>A</tt></li>
</ul>

<h2><a name="misc_fn"></a>3.36 Miscellaneous functions</h2>
<p class="small"><A HREF="#error">next</A> - <A HREF="#error">skip</A> -
<A HREF="#refer">up</A> - <A HREF="#top">start</A></p>
<p>The section includes some miscellaneous functions I have needed for my work. So 
far there are only the Helmert transforms.</p>
<h4>Helmert transforms</h4>
<p>This section refers to the Helmert transform used in some statistical 
packages for extracting contrasts. It is different from the Helmert transform 
used in geodesy. The version of the transform I am going to use is based on the 
<tt>n x n</tt> matrix</p>
<pre class="small">| 1/sqrt(1*2)                       | * | -1  1              |
|     1/sqrt(2*3)                   |   | -1 -1  2           |
|          1/sqrt(3*4)              |   | -1 -1 -1  3        |
|               ...                 |   |  ...               |
|                  1/(sqrt(n-1)*n)  |   | -1 -1  ...  -1 n-1 |
|                         1/sqrt(n) |   |  1  1  ...       1 | </pre>
<p>You can interpret multiplying a column vector by this matrix as follows. Form 
the <i>k</i>-th element in the resulting column vector by taking the <i>k+1</i>-th 
element and subtracting the average of the previous <i>k</i> elements. Then 
multiply by sqrt((<i>k+1</i>)/<i>k</i>) so we get an orthonormal transform. The 
last element is just the sum divided by sqrt(<i>n</i>). Usually one will omit 
the last element since we want just the contrasts.</p>
<p>Here are the functions. <tt>X</tt> and <tt>Y</tt> are ColumnVectors, <tt>A</tt> 
and <tt>B</tt> are matrices, <tt>b</tt> is a boolean, <tt>j</tt>, <tt>n</tt> are 
integers, <tt>x</tt> is a Real.</p>


<TABLE WIDTH=100% id="table1">
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A = Helmert(n, b);</tt><br>
<tt>A = Helmert(n);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Return the <i>n</i> x <i>n</i> 
Helmert matrix or the (<i>n</i>-1) x <i>n</i> version with the last row omitted.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>Y = Helmert(X,b);</tt><br>
<tt>Y = Helmert(X);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Multiply by the Helmert matrix.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>Y = Helmert(n,j,b);</tt><br>
<tt>Y = Helmert(n,j);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Return <i>j</i>-th column of Helmert 
matrix as a ColumnVector.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>X = Helmert_transpose(Y,b);</tt><br>
<tt>X = Helmert_transpose(Y);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Multiply by the transpose of the 
Helmert matrix.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>x = Helmert_transpose(Y,j,b);</tt><br>
<tt>x = Helmert_transpose(Y,j);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Multiply by the transpose of the 
Helmert matrix, return just the <i>j</i>-th element.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>B = Helmert(A,b);</tt><br>
<tt>B = Helmert(A);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Multiply a Matrix by the Helmert 
matrix.</TD>
</TR>
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <tt>A = Helmert_transpose(B,b);</tt><br>
<tt>A = Helmert_transpose(B);</tt></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> Multiply a Matrix by the transpose of 
the Helmert matrix.</TD>
</TR>
</TABLE>

<p>If <i>b</i> is true the full Helmert matrix is used, if it is false or 
omitted the <i>n</i>-th row (or <i>n</i>-column of the transpose) is omitted.</p>
<p>&nbsp;</p>

<H2><A NAME="error"></A>4. Error messages</H2>
<P CLASS="small"><A HREF="#design">next</A> - <A HREF="#design">skip</A> -
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
<P>Most error messages are self-explanatory. The message gives the size of the
matrices involved. Matrix types are referred to by the following codes: </P>
<PRE>   Matrix or vector                   Rect
   Symmetric matrix                   Sym
   Band matrix                        Band
   Symmetric band matrix              SmBnd
   Lower triangular matrix            LT
   Lower triangular band matrix       LwBnd
   Upper triangular matrix            UT
   Upper triangular band matrix       UpBnd
   Diagonal matrix                    Diag
   Crout matrix (LU matrix)           Crout
   Band LU matrix                     BndLU
</PRE>

<P>Other codes should not occur. </P>
<P>See the section on <A HREF="#except">exceptions</A> for more details on the
structure of the exception classes. </P>
<P>I have defined a class Tracer that is intended to help locate the place
where an error has occurred. At the beginning of a function I suggest you
include a statement like </P>
<PRE>   Tracer tr(&quot;name&quot;);
</PRE>

<P>where name is the name of the function. This name will be printed as part of
the error message, if an exception occurs in that function, or in a function
called from that function. You can change the name as you proceed through a
function with the ReName function </P>
<PRE>   tr.ReName(&quot;new name&quot;);
</PRE>

<P>if, for example, you want to track progress through the function. </P>
<H2><A NAME="design"></A>5. Notes on the design of the
library</H2>
<P CLASS="small"><A HREF="#sue">next</A> - <a href="#function">skip</a> -
<A HREF="#top">up</A> - <A HREF="#top">start</A></P>
<TABLE WIDTH="100%">
<TR>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A HREF="#sue">5.1
Safety, usability, efficiency </A><BR>
<A HREF="#mat_arr">5.2 Matrix vs array library </A><BR>
<A HREF="#question">5.3 Design questions </A><BR>
<A HREF="#stor">5.4 Data storage </A><BR>
<A HREF="#mem_man">5.5 Memory management - 1 </A><BR>
<A HREF="#mem_man2">5.6 Memory management - 2 </A><BR>
<A HREF="#evalx">5.7 Evaluation of expressions </A></TD>
<TD VALIGN="TOP" ALIGN="LEFT" WIDTH="50%"> <A
HREF="#explode">5.8 Explosion in the number of operations </A><BR>
<A HREF="#destr">5.9 Destruction of temporaries </A><BR>
<A HREF="#calc">5.10 A calculus of matrix types </A><BR>
<A HREF="#pointer">5.11 Pointer arithmetic </A><BR>
<A HREF="#err_hand">5.12 Error handling </A><BR>
<A HREF="#sparse">5.13 Sparse matrices </A><BR>
<A HREF="#comp_mat">5.14 Complex matrices </A></TD>
</TR>
</TABLE>
<P>I describe some of the ideas behind this package, some of the decisions that
I needed to make and give some details about the way it works. You don't need
to read this part of the documentation in order to use the package. </P>
<P>It isn't obvious what is the best way of going about structuring a matrix
package. I don't think you can figure this out with <I>thought</I> experiments.
Different people have to try out different approaches. And someone else may
have to figure out which is best. Or, more likely, the ultimate packages will
lift some ideas from each of a variety of trial packages. So, I don't claim my
package is an <I>ultimate</I> package, but simply a trial of a number of ideas.
The following pages give some background on these ideas. </P>
<H2><A NAME="sue"></A>5.1 Safety, usability,
efficiency</H2>
<P CLASS="small"><A HREF="#mat_arr">next</A> - <A HREF="#mat_arr">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
<H3>Some general comments</H3>
<P>A library like <I>newmat</I> needs to balance <I>safety</I>,
<I>usability</I> and <I>efficiency</I>. </P>
<P>By <B>safety</B>, I mean getting the right answer, and not causing crashes
or damage to the computer system. </P>
<P>By <B>usability</B>, I mean being easy to learn and use, including not being
too complicated, being intuitive, saving the users' time, being nice to use. 
</P>
<P><B>Efficiency</B> means minimising the use of computer memory and time. </P>
<P>In the early days of computers the emphasis was on efficiency. But computer
power gets cheaper and cheaper, halving in price every 18 months. On the other
hand the unaided human brain is probably not a lot better than it was 100,000
years ago! So we should expect the balance to shift to put more emphasis on
safety and usability and a little less on efficiency. So I don't mind if my
programs are a little less efficient than programs written in pure C (or
Fortran) if I gain substantially in safety and usability. But I would mind if
they were a lot less efficient. </P>
<H3>Type of use</H3>
<P>Second reason for putting extra emphasis on safety and usability is the way
I and, I suspect, most other users actually use <I>newmat</I>. Most completed
programs are used only a few times. Some result is required for a client, paper
or thesis. The program is developed and tested, the result is obtained, and the
program archived. Of course bits of the program will be recycled for the next
project. But it may be less usual for the same program to be run over and over
again. So the cost, computer time + people time, is in the development time and
often, much less in the actual time to run the final program. So good use of
people time, especially during development is really important. This means you
need highly usable libraries. </P>
<P>So if you are dealing with matrices, you want the good interface that I have
tried to provide in <I>newmat</I>, and, of course, reliable methods underneath
it. </P>
<P>Of course, efficiency is still important. We often want to run the biggest
problem our computer will handle and often a little bigger. The C++ language
almost lets us have both worlds. We can define a reasonably good interface, and
get good efficiency in the use of the computer. </P>
<H3>Levels of access</H3>
<P>We can imagine the <I>black box</I> model of a <I>newmat</I>. Suppose the
inside is hidden but can be accessed by the methods described in the
<A HREF="#refer">reference</A> section. Then the interface is reasonably
consistent and intuitive. Matrices can be accessed and manipulated in much the
same way as doubles or ints in regular C. All accesses are checked. It is most
unlikely that an incorrect index will crash the system. In general, users do
not need to use pointers, so one shouldn't get pointers pointing into space.
And, hopefully, you will get simpler code with less errors. </P>
<P>There are some exceptions to this. In particular, the <A
HREF="#elements">C-like subscripts</A> are not checked for validity. They give
faster access but with a lower level of safety. </P>
<P>Then there is the <a href="#scalar1">data()</a> function which takes you to
the data array within a matrix. This takes you right inside the <I>black
box</I>. But this is what you have to use if you are writing, for example, a
new matrix factorisation, and require fast access to the data array. I have
tried to write code to simplify access to the interior of a rectangular matrix,
see file newmatrm.cpp, but I don't regard this as very successful, as yet, and
have not included it in the documentation. Ideally we should have improved
versions of this code for each of the major types of matrix. But, in reality,
most of my matrix factorisations are written in what is basically the C
language with very little C++. </P>
<P>So our <I>box</I> is not very <I>black</I>. You have a choice of how far you
penetrate. On the outside you have a good level of safety, but in some cases
efficiency is compromised a little. If you penetrate inside the <I>box</I>
safety is reduced but you can get better efficiency. </P>
<H3>Some performance data</H3>
<P>This section looks at the performance on <I>newmat</I> for simple sums,
comparing it with C code and with a simple array program. 
(This is now very old data).</P>
<P>The following table lists the time (in seconds) for carrying out the
operations <TT>X=A+B;</TT>, <TT>X=A+B+C;</TT>, <TT>X=A+B+C+D;</TT>, <TT>X=A+B+C+D+E;</TT> where
<TT>X,A,B,C,D,E</TT> are of type ColumnVector, with a variety of programs. I am
using Microsoft VC++, version 6 in  console mode under windows 2000 on a PC with 
a 1 ghz Pentium III and 512 mbytes of memory. </P>
<pre>    length    iters. newmat      C C-res.  subs.  array
<b>X = A + B</b>
         2   5000000   27.8    0.3    8.8    1.9    9.5 
        20    500000    3.0    0.3    1.1    1.9    1.2 
       200     50000    0.5    0.3    0.4    1.9    0.3 
      2000      5000    0.4    0.3    0.4    2.0    1.0 
     20000       500    4.5    4.5    4.5    6.7    4.4 
    200000        50    5.2    4.7    5.5    5.8    5.2 

<b>X = A + B + C</b>
         2   5000000   36.6    0.4    8.9    2.5   12.2 
        20    500000    4.0    0.4    1.2    2.5    1.6 
       200     50000    0.8    0.3    0.5    2.5    0.5 
      2000      5000    3.6    4.4    4.6    9.0    4.4 
     20000       500    6.8    5.4    5.4    9.6    6.8 
    200000        50    8.6    6.0    6.7    7.1    8.6 

<b>X = A + B + C + D</b>
         2   5000000   44.0    0.7    9.3    3.1   14.6 
        20    500000    4.9    0.6    1.5    3.1    1.9 
       200     50000    1.0    0.6    0.8    3.2    0.8 
      2000      5000    5.6    6.6    6.8   11.5    5.9 
     20000       500    9.0    6.7    6.8   11.0    8.5 
    200000        50   11.9    7.1    7.9    9.5   12.0 

<b>X = A + B + C + D + E</b>
         2   5000000   50.6    1.0    9.5    3.8   17.1 
        20    500000    5.7    0.8    1.7    3.9    2.4 
       200     50000    1.3    0.9    1.0    3.9    1.0 
      2000      5000    7.0    8.3    8.2   13.8    7.1 
     20000       500   11.5    8.1    8.4   13.2   11.0 
    200000        50   15.2    8.7    9.5   12.4   15.4 
</pre>
<P>I have <a target="_blank" href="add_time.png">graphed</a> the results 
and included rather more array lengths.</P>

<P>The first column gives the lengths of the arrays, the second the number of
iterations and the remaining columns the total time required in seconds. If the
only thing that consumed time was the double precision addition then the
numbers within each block of the table would be the same. The summation is 
repeated 5 times within each loop, for example:</P>
<pre>   for (i=1; i&lt;=m; ++i)
   {
      X1 = A1+B1+C1; X2 = A2+B2+C2; X3 = A3+B3+C3;
      X4 = A4+B4+C4; X5 = A5+B5+C5;
   }</pre>
<P>The column labelled <I>newmat</I> is using the standard <I>newmat</I> add. The column labelled <I>C</I> uses the usual C method: <TT>while
(j1--) *x1++ = *a1++ + *b1++;</TT> . The following column also includes an
<TT>X.resize()</TT> in the outer loop to correspond to the reassignment of
memory that <I>newmat</I> would do. In the next column the calculation is using 
the usual C style <I>for</I> loop
and accessing the elements using <I>newmat</I> subscripts such as
<TT>A(i)</TT>. The final column is the time taken by a
simple array package. This uses an alternative method for avoiding temporaries 
and unnecessary copies that does not involve runtime tests. It does its sums in blocks of 4 and copies in blocks of
8 in the same way that <I>newmat</I> does. </P>
<P>Here are my conclusions. </P>
<UL>
<LI><I>Newmat</I> does very badly for length 2 and doesn't do  well for
length 20. There is  a lot of code in <I>newmat</I> for
determining which sum algorithm to use and it is not surprising that this 
impacts on performance for small lengths.
However the <I>array</I> program is also having difficulty with length 2 so it
is unlikely that the problem could be completely eliminated. </LI>
<LI>For arrays of length 2000 or longer <I>newmat</I> is doing about as well as
C and slightly better than C with resize in the <TT>X=A+B</TT> table. For the
other two tables it tends to be slower, but not dramatically so. </LI>
<LI>It is really important for fast processing with the Pentium III to stay 
within the Pentium cache.</LI>
<LI>Addition using the <I>newmat</I> subscripts, while considerably slower than
the others, is still surprisingly good for the longer arrays. </LI>
<LI>The <I>array</I> program and <I>newmat</I> are similar for
lengths 2000 or higher (the longer times for the array program for the longest 
arrays shown on the graph are probably a quirk of the timing program). </LI>
</UL>
<P>In summary: for the situation considered here, <I>newmat</I> is doing very
well for large ColumnVectors, even for sums with several terms, but not so well
for shorter ColumnVectors. </P>
<H2><A NAME="mat_arr"></A>5.2 Matrix vs array
library</H2>
<P CLASS="small"><A HREF="#question">next</A> - <A HREF="#question">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>The <I>newmat</I> library is for the manipulation of matrices, including the
standard operations such as multiplication as understood by numerical analysts,
engineers and mathematicians. </P>
<P>A matrix is a two dimensional array of numbers. However, very special
operations such as matrix multiplication are defined specifically for matrices.
This means that a <I>matrix</I> library, as I understand the term, is different
from a general <I>array</I> library. Here are some contrasting properties.</P>
<TABLE WIDTH="100%" BORDER="1">
<TR>
<TH VALIGN="TOP" WIDTH="20%" ALIGN="LEFT">Feature</TH>
<TH VALIGN="TOP" WIDTH="40%" ALIGN="LEFT">Matrix
library</TH>
<TH VALIGN="TOP" WIDTH="40%" ALIGN="LEFT">Array
library</TH>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Expressions</TD>
<TD VALIGN="TOP" WIDTH="40%">Matrix expressions<TT>;</TT>
<TT>*</TT> means matrix multiply; inverse function</TD>
<TD VALIGN="TOP" WIDTH="40%">Arithmetic operations, if
supported, mean elementwise combination of arrays</TD>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Element access</TD>
<TD VALIGN="TOP" WIDTH="40%">Access to the elements of a
matrix</TD>
<TD VALIGN="TOP" WIDTH="40%">High speed access to elements
directly and perhaps with iterators</TD>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Elementary functions</TD>
<TD VALIGN="TOP" WIDTH="40%">For example: determinant,
trace</TD>
<TD VALIGN="TOP" WIDTH="40%">Matrix multiplication as a
function</TD>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Advanced functions</TD>
<TD VALIGN="TOP" WIDTH="40%">For example: eigenvalue
analysis</TD>
<TD VALIGN="TOP" WIDTH="40%">&nbsp;</TD>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Element types</TD>
<TD VALIGN="TOP" WIDTH="40%">Real and possibly
complex</TD>
<TD VALIGN="TOP" WIDTH="40%">Wide range - real, integer, string
etc</TD>
</TR>
<TR>
<TD VALIGN="TOP" WIDTH="20%">Types</TD>
<TD VALIGN="TOP" WIDTH="40%">Rectangular, symmetric, diagonal,
etc</TD>
<TD VALIGN="TOP" WIDTH="40%">One, two and three dimensional
arrays, at least</TD>
</TR>
</TABLE>
<P>Both types of library need to support access to sub-matrices or sub-arrays,
have good efficiency and storage management, and graceful exit for errors. In
both cases, we probably need two versions, one optimised for large matrices or
arrays and one for small matrices or arrays.</P>
<P>It may be possible to amalgamate the two sets of requirements to some
extent. However <I>newmat</I> is definitely oriented towards the matrix library
set.</P>
<H2><A NAME="question"></A>5.3 Design questions</H2>
<P CLASS="small"><A HREF="#stor">next</A> - <A HREF="#stor">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>Even within the bounds set by the requirements of a matrix library there is
a substantial opportunity for variation between what different matrix packages
might provide. It is not possible to build a matrix package that will meet
everyone's requirements. In many cases if you put in one facility, you impose
overheads on everyone using the package. This both in storage required for the
program and in efficiency. Likewise a package that is optimised towards
handling large matrices is likely to become less efficient for very small
matrices where the administration time for the matrix may become significant
compared with the time to carry out the operations. It is better to provide a
variety of packages (hopefully compatible) so that most users can find one that
meets their requirements. This package is intended to be one of these packages;
but not all of them. </P>
<P>Since my background is in statistical methods, this package is oriented
towards the kinds things you need for statistical analyses. </P>
<P>Now looking at some specific questions. </P>
<H3>What size of matrices?</H3>
<P>A matrix library may target small matrices (say 3 x 3), or medium sized
matrices, or very large matrices. </P>
<P>A library targeting very small matrices will seek to minimise
administration. A library for medium sized or very large matrices can spend
more time on administration in order to conserve space or optimise the
evaluation of expressions. A library for very large matrices will need to pay
special attention to storage and numerical properties. This library is designed
for medium sized matrices. This means it is worth introducing some
optimisations, but I don't have to worry about setting up some form of virtual
memory. </P>
<H3>Which matrix types?</H3>
<P>As well as the usual rectangular matrices, matrices occurring repeatedly in
numerical calculations are upper and lower triangular matrices, symmetric
matrices and diagonal matrices. This is particularly the case in calculations
involving least squares and eigenvalue calculations. So as a first stage these
were the types I decided to include. </P>
<P>It is also necessary to have types row vector and column vector. In a
<I>matrix</I> package, in contrast to an <I>array</I> package, it is necessary
to have both these types since they behave differently in matrix expressions.
The vector types can be derived for the rectangular matrix type, so having them
does not greatly increase the complexity of the package. </P>
<P>The problem with having several matrix types is the number of versions of
the binary operators one needs. If one has 5 distinct matrix types then a
simple library will need 25 versions of each of the binary operators. In fact,
we can evade this problem, but at the cost of some complexity. </P>
<H3>What element types?</H3>
<P>Ideally we would allow element types double, float, complex and int, at
least. It might be reasonably easy, using templates or equivalent, to provide a
library which could handle a variety of element types. However, as soon as one
starts implementing the binary operators between matrices with different
element types, again one gets an explosion in the number of operations one
needs to consider. At the present time the compilers I deal with are not up to
handling this problem with templates. (Of course, when I started writing
<I>newmat</I> there were no templates). But even when the compilers do meet the
specifications of the draft standard, writing a matrix package that allows for
a variety of element types using the template mechanism is going to be very
difficult. I am inclined to use templates in an <I>array</I> library but not in
a <I>matrix</I> library. </P>
<P>Hence I decided to implement only one element type. But the user can decide
whether this is float or double. The package assumes elements are of type Real.
The user typedefs Real to float or double. </P>
<P>It might also be worth including symmetric and triangular matrices with
extra precision elements (double or long double) to be used for storage only
and with a minimum of operations defined. These would be used for accumulating
the results of sums of squares and product matrices or multi-stage QR 
decompositions. </P>
<H3>Allow matrix expressions</H3>
<P>I want to be able to write matrix expressions the way I would on paper. So
if I want to multiply two matrices and then add the transpose of a third one I
can write something like <TT>X = A * B + C.t();</TT>. I want this expression to
be evaluated with close to the same efficiency as a hand-coded version. This is
not so much of a problem with expressions including a multiply since the
multiply will dominate the time. However, it is not so easy to achieve with
expressions with just <TT>+</TT> and <TT>-</TT>. </P>
<P>A second requirement is that temporary matrices generated during the
evaluation of an expression are destroyed as quickly as possible. </P>
<P>A desirable feature is that a certain amount of <I>intelligence</I> be
displayed in the evaluation of an expression. For example, in the expression
<TT>X = A.i() * B;</TT> where <TT>i()</TT> denotes inverse, it would be
desirable if the inverse wasn't explicitly calculated. </P>
<H3>Naming convention</H3>
<P>How are classes and public member functions to be named? As a general rule I
have spelt identifiers out in full with individual words being capitalised. For
example <I>UpperTriangularMatrix</I>. If you don't like this you can #define or
typedef shorter names. This convention means you can select an abbreviation
scheme that makes sense to you. </P>
<P>Exceptions to the general rule are the functions for transpose and inverse.
To make matrix expressions more like the corresponding mathematical formulae, I
have used the single letter abbreviations, <TT>t()</TT> and <TT>i()</TT>. </P>
<P>I am now switching to using lowercase for  functions with individual 
words separated by &quot;<tt>_</tt>&quot;. This is following the convention in 
the standard library and I think it looks neater. Class names will remain with 
individual words being capitalised.</P>
<H3>Row and column index ranges</H3>
<P>In mathematical work matrix subscripts usually start at one. In C, array
subscripts start at zero. In Fortran, they start at one. Possibilities for this
package were to make them start at 0 or 1 or be arbitrary. </P>
<P>Alternatively one could specify an <I>index set</I> for indexing the rows
and columns of a matrix. One would be able to add or multiply matrices only if
the appropriate row and column index sets were identical. </P>
<P>In fact, I adopted the simpler convention of making the rows and columns of
a matrix be indexed by an integer starting at one, following the traditional
convention. In an earlier version of the package I had them starting at zero,
but even I was getting mixed up when trying to use this earlier package. So I
reverted to the more usual notation and started at 1. </P>
<H3>Element access - method and checking</H3>
<P>We want to be able to use the notation <TT>A(i,j)</TT> to specify the
<TT>(i,j)</TT>-th element of a matrix. This is the way mathematicians expect to
address the elements of matrices. I consider the notation <TT>A[i][j]</TT>
totally alien. However I include this as an option to help people converting
from C. </P>
<P>There are two ways of working out the address of <TT>A(i,j)</TT>. One is
using a <I>dope</I> vector which contains the first address of each row.
Alternatively you can calculate the address using the formula appropriate for
the structure of <TT>A</TT>. I use this second approach. It is probably slower,
but saves worrying about an extra bit of storage. </P>
<P>The other question is whether to check for <TT>i</TT> and <TT>j</TT> being
in range. I do carry out this check following years of experience with both
systems that do and systems that don't do this check. I would hope that the
routines I supply with this package will reduce your need to access elements of
matrices so speed of access is not a high priority. </P>
<H3>Use iterators</H3>
<P>Iterators are an alternative way of providing fast access to the elements of
an array or matrix when they are to be accessed sequentially. They need to be
customised for each type of matrix. I have not implemented iterators in this
package, although some iterator like functions are used internally for some row
and column functions. </P>
<H2><A NAME="stor"></A>5.4 Data storage</H2>
<P CLASS="small"><A HREF="#mem_man">next</A> - <A HREF="#mem_man">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<H3>The stack and heap</H3>
<p>To understand how <i>newmat</i> stores matrices you need to know a little bit 
about the <i>heap</i> and <i>stack</i>.</p>
<p>The data values of variables or objects in a C++ program are stored in either 
of two sections of memory called the <i>stack</i> and the <i>heap</i>. Sometimes 
there is more than one <i>heap</i> to cater for different sized variables.</p>
<p>If you declare an <i>automatic</i> variable</p>
<pre>   int x;</pre>
<p>then the value of <i>x</i> is stored on the <i>stack</i>. As you declare more 
variables the stack gets bigger. When you exit a block (i.e a section of code 
delimited by curly brackets <i>{...}</i>) the memory used by the automatic 
variables declared in the block is released and the <i>stack</i> shrinks.</p>
<p>When you declare a variable with <i>new</i>, for example,</p>
<pre>   int* y = new int;</pre>
<p>the pointer <i>y</i> is stored on the <i>stack</i> but the value it is 
pointing to is stored on the <i>heap</i>. Memory on the <i>heap</i> is not 
released until the program explicitly does this with a <i>delete</i> statement</p>
<pre>   delete y;</pre>
<p>or the program exits.</p>
<p>On the <i>stack</i>, variables and objects are is always added to the end of 
the <i>stack</i> and are removed in the reverse order to that in which they are 
added - that is the last on will be the first off. This is not the case with the <i>
heap</i>, where the variables and objects can be removed in any order. So one 
can get alternating pieces of used and unused memory. When a new variable or 
object is declared on the <i>heap</i> the system needs to search for piece of 
unused memory large enough to hold it. This means that storing on the <i>heap</i> 
will usually be a slower process than storing on the <i>stack</i>. There is also 
likely to be waste space on the <i>heap</i> because of gaps between the used 
blocks of memory that are too small for the next object you want to store on the
<i>heap</i>. There is also the possibility of wasting space if you forget to 
remove a variable or object on the <i>heap</i> even though you have finished 
using it. However, the <i>stack</i> is usually limited to holding small objects 
with size known at compile time. Large objects, objects whose size you don't 
know at compile time, and objects that you want to persist after the end of the 
block need to be stored on the <i>heap</i>.</p>
<p>In C++, the <i>constructor</i>/<i>destructor</i> system enables one to build 
complicated objects such as matrices that behave as automatic variables stored 
on the <i>stack</i>, so the programmer doesn't have to worry about deleting them 
at the end of the block, but which really utilise the <i>heap</i> for storing 
their data.</p>
<H3>Structure of matrix objects</H3>
<P>Each matrix object contains the basic information such as the number of rows
and columns, the amount of memory used, a status variable and a pointer to the data array which is on
the heap. So if you declare a matrix</P>
<pre>   Matrix A(1000,1000);</pre>
<p>there is an small amount of memory used on the stack for storing the numbers 
of rows and columns, the amount of&nbsp; memory used, the status variable and 
the pointer together with 1,000,000 <i>Real</i> locations stored on the heap. 
When you exit the block in which <i>A</i> is declared, the heap memory used by
<i>A</i> is automatically returned to the system, as well as the memory used on 
the stack.</p>
<p>Of course, if you use new to declare a matrix</p>
<pre>   Matrix* B = new Matrix(1000,1000);</pre>
<p>both the information about the size and the actual data are stored on heap 
and not deleted until the program exits or you do an explicit delete:</p>
<pre>   delete B;</pre>
<p>If you carry out an assignment with <tt>=</tt> or <tt>&lt;&lt;</tt> or do a 
<tt>resize()</tt> the data array currently associated with a matrix is destroyed and 
a new array generated. For example</p>
<pre>   Matrix A(1000,1000);
   Matrix B(50, 50);
   ... put some values in B
   A = B;</pre>
<p>At the last step the heap memory associated with <i>A</i> is returned to the 
system and a new block of heap memory is assigned to contain the new values. 
This happens even if there is no change in the amount of memory required. </p>
<H3>One block or several</H3>
<P>The elements of the matrix are stored as a single array. Alternatives would
have been to store each row as a separate array or a set of adjacent rows as a
separate array. The present solution simplifies the program but limits the size
of matrices in 16 bit PCs that have a 64k byte limit on the size of arrays (I
don't use the <TT>huge</TT> keyword). The large arrays may also cause problems
for memory management in smaller machines. [The 16 bit PC problem has largely
gone away but it was a problem when much of <I>newmat</I> was written. Now,
occasionally I run into the 32 bit PC problem.]</P>
<H3>By row or by column or other</H3>
<P>In Fortran two dimensional arrays are stored by column. In most other
systems they are stored by row. I have followed this later convention. This
makes it easier to interface with other packages written in C but harder to
interface with those written in Fortran. This may have been a wrong decision.
Most work on the efficient manipulation of large matrices is being done in
Fortran. It would have been easier to use this work if I had adopted the
Fortran convention. </P>
<P>An alternative would be to store the elements by mid-sized rectangular
blocks. This might impose less strain on memory management when one needs to
access both rows and columns. </P>
<H3>Storage of symmetric matrices</H3>
<P>Symmetric matrices are stored as lower triangular matrices. The decision was
pretty arbitrary, but it does slightly simplify the Cholesky decomposition
program. </P>
<H2><A NAME="mem_man"></A>5.5 Memory management -
reference counting or status variable?</H2>
<P CLASS="small"><A HREF="#mem_man2">next</A> - <A HREF="#mem_man2">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>Consider the instruction </P>
<PRE>   X = A + B + C;
</PRE>

<P>To evaluate this a simple program will add <TT>A</TT> to <TT>B</TT> putting
the total in a temporary <TT>T1</TT>. Then it will add <TT>T1</TT> to
<TT>C</TT> creating another temporary <TT>T2</TT> which will be copied into
<TT>X</TT>. <TT>T1</TT> and <TT>T2</TT> will sit around till the end of the
execution of the statement and perhaps of the block. It would be faster if the
program recognised that <TT>T1</TT> was temporary and stored the sum of
<TT>T1</TT> and <TT>C</TT> back into <TT>T1</TT> instead of creating
<TT>T2</TT> and then avoided the final copy by just assigning the contents of
<TT>T1</TT> to <TT>X</TT> rather than copying. In this case there will be no
temporaries requiring deletion. (More precisely there will be a header to be
deleted but no contents). </P>
<P>For an instruction like </P>
<PRE>   X = (A * B) + (C * D);
</PRE>

<P>we can't easily avoid one temporary being left over, so we would like this
temporary deleted as quickly as possible. </P>
<P>I provide the functionality for doing all this by attaching a status
variable to each matrix. This indicates if the matrix is temporary so that its
memory is available for recycling or deleting. Any matrix operation checks the
status variables of the matrices it is working with and recycles or deletes any
temporary memory. </P>
<P>An alternative or additional approach would be to use <I>reference counting
and delayed copying</I> - also known as <I>copy on write</I>. If a program
requests a matrix to be copied, the copy is delayed until an instruction is
executed which modifies the memory of either the original matrix or the copy.
If the original matrix is deleted before either matrix is modified, in effect,
the values of the original matrix are transferred to the copy without any
actual copying taking place. This solves the difficult problem of returning an
object from a function without copying and saves the unnecessary copying in the
previous examples. </P>
<P>There are downsides to the delayed copying approach. Typically, for delayed
copying one uses a structure like the following: </P>
<PRE>   Matrix
     |
     +------&gt; Array Object
     |          |
     |          +------&gt; Data array
     |          |
     |          +------- Counter
     |
     +------ Dimension information

</PRE>

<P>where the arrows denote a pointer to a data structure. If one wants to
access the <I>Data array</I> one will need to track through two pointers. If
one is going to write, one will have to check whether one needs to copy first.
This is not important when one is going to access the whole array, say, for a
add operation. But if one wants to access just a single element, then it
imposes a significant additional overhead on that operation. Any subscript
operation would need to check whether an update was required - even read since
it is hard for the compiler to tell whether a subscript access is a read or
write. </P>
<P>Some matrix libraries don't bother to do this. So if you write <TT>A =
B;</TT> and then modify an element of one of <TT>A</TT> or <TT>B</TT>, then the
same element of the other is also modified. I don't think this is acceptable
behaviour. </P>
<P>Delayed copy does not provide the additional functionality of my approach
but I suppose it would be possible to have both delayed copy and tagging
temporaries. </P>
<P>My approach does not automatically avoid all copying. In particular, you
need use a special technique to return a matrix from a function without
copying. </P>
<H2><A NAME="mem_man2"></A>5.6 Memory management -
accessing contiguous locations</H2>
<P CLASS="small"><A HREF="#evalx">next</A> - <A HREF="#evalx">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>Modern computers work faster if one accesses memory by running through
contiguous locations rather than by jumping around all over the place. Newmat
stores matrices <A HREF="#stor">by rows</A> so that algorithms that access
memory by running along rows will tend to work faster than one that runs down
columns. A number of the algorithms used in <i>Newmat</i> were developed before this
was an issue and so are not as efficient as possible. </P>
<P>I have gradually upgrading the algorithms to access memory by rows. The
following table shows the current status of this process.</P>
<TABLE WIDTH="100%" BORDER="1">
<TR>
<TH ALIGN="LEFT">Function</TH>
<TH ALIGN="LEFT">Contiguous memory access</TH>
<TH ALIGN="LEFT">Comment</TH>
</TR>
<TR>
<TD>Add, subtract</TD>
<TD>Yes</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Multiply</TD>
<TD>Yes</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Concatenate</TD>
<TD>Yes</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Transpose</TD>
<TD>No</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Invert and solve</TD>
<TD>Yes</TD>
<TD>Mostly</TD>
</TR>
<TR>
<TD>Cholesky</TD>
<TD>Yes</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>QRZ, QRZT</TD>
<TD>Yes</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>updateQRZ, updateQRZT</TD>
<TD>Partially</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>SVD</TD>
<TD>No</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Jacobi</TD>
<TD>No</TD>
<TD>Not an issue; used only for smaller matrices</TD>
</TR>
<TR>
<TD>Eigenvalues</TD>
<TD>No</TD>
<TD>&nbsp;</TD>
</TR>
<TR>
<TD>Sort</TD>
<TD>Yes</TD>
<TD>Quick-sort is naturally good</TD>
</TR>
<TR>
<TD>FFT</TD>
<TD>?</TD>
<TD>Could be improved?</TD>
</TR>
</TABLE>
<p>This is now all rather out of date. With Pentiums, at least, the important 
requirement for speed seems to be to minimise transfers between the RAM memory 
and the on-chip memory. There isn't much you can do about add and subtract, but 
there lots of possibilities for some of the other operations.</p>
<H2><A NAME="evalx"></A>5.7 Evaluation of expressions -
lazy evaluation</H2>
<P CLASS="small"><A HREF="#explode">next</A> - <A HREF="#explode">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
<P>Consider the instruction </P>
<PRE>   X = B - X;
</PRE>

<P>A simple program will subtract <TT>X</TT> from <TT>B</TT>, store the result
in a temporary <TT>T1</TT> and copy <TT>T1</TT> into <TT>X</TT>. It would be
faster if the program recognised that the result could be stored directly into
<TT>X</TT>. This would happen automatically if the program could look at the
instruction first and mark <TT>X</TT> as temporary. </P>
<P>C programmers would expect to avoid the same problem with </P>
<PRE>   X = X - B;
</PRE>

<P>by using an operator <TT>-=</TT> </P>
<PRE>   X -= B;
</PRE>

<P>However this is an unnatural notation for non C users and it may be nicer to
write <TT>X = X - B</TT>; and know that the program will carry out the
simplification. </P>
<P>Another example where this intelligent analysis of an instruction is helpful
is in </P>
<PRE>   X = A.i() * B;
</PRE>

<P>where <TT>i()</TT> denotes inverse. Numerical analysts know it is
inefficient to evaluate this expression by carrying out the inverse operation
and then the multiply. Yet it is a convenient way of writing the instruction.
It would be helpful if the program recognised this expression and carried out
the more appropriate approach. </P>
<P>I regard this interpretation of <TT>A.i() * B</TT> as just providing a
convenient notation. The objective is not primarily to correct the errors of
people who are unaware of the inefficiency of <TT>A.i() * B</TT> if interpreted
literally. </P>
<P>There is a third reason for the two-stage evaluation of expressions and this
is probably the most important one. In C++ it is quite hard to return an
expression from a function such as (<TT>*</TT>, <TT>+</TT> etc) without a copy.
This is particularly the case when an assignment (<TT>=</TT>) is involved. The
mechanism described here provides one way for avoiding this in matrix
expressions. </P>
<P>The C++ standard (section 12.8/15) allows the compiler to optimise away the
copy when returning an object from a function (but there will still be one copy
is an <span lang="en-nz">assignment </span>(=) is involved). This means special handling of returns from a
function is less important when a modern optimising compiler is being
used.&nbsp; </P>
<P>To carry out this <I>intelligent</I> analysis of an instruction matrix
expressions are evaluated in two stages. In the the first stage a tree
representation of the expression is formed. For example <TT>(A+B)*C</TT> is
represented by a tree </P>
<PRE><TT>
       *
      / \
     +   C
    / \
   A   B
</TT></PRE>

<P>Rather than adding <TT>A</TT> and <TT>B</TT> the <TT>+</TT> operator yields
an object of a class <I>AddedMatrix</I> which is just a pair of pointers to
<TT>A</TT> and <TT>B</TT>. Then the <TT>*</TT> operator yields a
<I>MultipliedMatrix</I> which is a pair of pointers to the <I>AddedMatrix</I>
and <TT>C</TT>. The tree is examined for any simplifications and then evaluated
recursively. </P>
<P>Further possibilities not yet included are to recognise <TT>A.t()*A</TT> and
<TT>A.t()+A</TT> as symmetric or to improve the efficiency of evaluation of
expressions like <TT>A+B+C</TT>, <TT>A*B*C</TT>, <TT>A*B.t()</TT> (<TT>t()</TT>
denotes transpose). </P>
<P>One of the disadvantages of the two-stage approach is that the types of
matrix expressions are determined at run-time. So the compiler will not detect
errors of the type </P>
<PRE>   Matrix M;
   DiagonalMatrix D;
   ....;
   D = M;
</PRE>

<P>We don't allow conversions using <TT>=</TT> when information would be lost.
Such errors will be detected when the statement is executed. </P>
<H2><A NAME="explode"></A>5.8 How to overcome an
explosion in number of operations</H2>
<P CLASS="small"><A HREF="#destr">next</A> - <A HREF="#destr">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>The package attempts to solve the problem of the large number of versions of
the binary operations required when one has a variety of types. </P>
<P>With <I>n</I> types of matrices the binary operations will each require
<I>n</I>-squared separate algorithms. Some reduction in the number may be
possible by carrying out conversions. However, the situation rapidly becomes
impossible with more than 4 or 5 types. Doug Lea told me that it was possible
to avoid this problem. I don't know what his solution is. Here's mine. </P>
<P>Each matrix type includes routines for extracting individual rows or
columns. I assume a row or column consists of a sequence of zeros, a sequence
of stored values and then another sequence of zeros. Only a single algorithm is
then required for each binary operation. The rows can be located very quickly
since most of the matrices are stored row by row. Columns must be copied and so
the access is somewhat slower. As far as possible my algorithms access the
matrices by row. </P>
<P>There is another approach. Each of the matrix types defined in this package
can be set up so both rows and columns have their elements at equal intervals
provided we are prepared to store the rows and columns in up to three chunks.
With such an approach one could write a single &quot;generic&quot; algorithm
for each of multiply and add. This would be a reasonable alternative to my
approach. </P>
<P>I provide several algorithms for operations like + . If one is adding two
matrices of the same type then there is no need to access the individual rows
or columns and a faster general algorithm is appropriate. </P>
<P>Generally the method works well. However symmetric matrices are not always
handled very efficiently (yet) since complete rows are not stored explicitly. 
</P>
<P>The original version of the package did not use this access by row or column
method and provided the multitude of algorithms for the combination of
different matrix types. The code file length turned out to be just a little
longer than the present one when providing the same facilities with 5 distinct
types of matrices. It would have been very difficult to increase the number of
matrix types in the original version. Apparently 4 to 5 types is about the
break even point for switching to the approach adopted in the present package. 
</P>
<P>However it must also be admitted that there is a substantial overhead in the
approach adopted in the present package for small matrices. The test program
developed for the original version of the package takes 30 to 50% longer to run
with the current version (though there may be some other reasons for this).
This is for matrices in the range 6x6 to 10x10. </P>
<P>To try to improve the situation a little I do provide an ordinary matrix
multiplication routine for the case when all the matrices involved are
rectangular. </P>
<H2><A NAME="destr"></A>5.9 Destruction of
temporaries</H2>
<P CLASS="small"><A HREF="#calc">next</A> - <A HREF="#calc">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>Versions before version 5 of newmat did not work correctly with Gnu C++
(version 5 or earlier). This was because the tree structure used to represent a
matrix expression was set up on the stack.&nbsp; Early versions of Gnu C++ destroyed temporary structures as soon as the
function that accesses them finished.  To overcome this problem, there was an 
option to store the temporaries forming the tree structure on the heap (created 
with new) and to delete them explicitly. Now that the C++ standards committee 
has said that temporary structures should not be destroyed before a statement 
finishes, I have deleted this option.</P>
<H2><A NAME="calc"></A>5.10 A calculus of matrix
types</H2>
<P CLASS="small"><A HREF="#pointer">next</A> - <A HREF="#pointer">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
<P>The program needs to be able to work out the class of the result of a matrix
expression. This is to check that a conversion is legal or to determine the
class of an intermediate result. To assist with this, a class MatrixType is
defined. Operators <TT>+</TT>, <TT>-</TT>, <TT>*</TT>, <TT>&gt;=</TT> are
defined to calculate the types of the results of expressions or to check that
conversions are legal. </P>
<P>Early versions of <I>newmat</I> stored the types of the results of
operations in a table. So, for example, if you multiplied an
UpperTriangularMatrix by a LowerTriangularMatrix, <I>newmat</I> would look up
the table and see that the result was of type Matrix. With this approach the
<A HREF="#explode">exploding</A> number of operations problem recurred although
not as seriously as when code had to be written for each pair of types. But
there was always the suspicion that somewhere, there was an error in one of
those 9x9 tables, that would be very hard to find. And the problem would get
worse as additional matrix types or operators were included. </P>
<P>The present version of <I>newmat</I> solves the problem by assigning
<I>attributes</I> such as <I>diagonal</I> or <I>band</I> or <I>upper
triangular</I> to each matrix type. Which attributes a matrix type has, is
stored as bits in an integer. As an example, the DiagonalMatrix type has the
bits corresponding to <I>diagonal</I>, <I>symmetric</I> and <I>band</I> equal
to 1. By looking at the attributes of each of the operands of a binary
operator, the program can work out the attributes of the result of the
operation with simple bitwise operations. Hence it can deduce an appropriate
type. The <I>symmetric</I> attribute is a minor problem because
<I>symmetric</I> * <I>symmetric</I> does not yield <I>symmetric</I> unless both
operands are <I>diagonal</I>. But otherwise very simple code can be used to
deduce the attributes of the result of a binary operation. </P>
<P>Tables of the types resulting from the binary operators are output at the
beginning of the <A HREF="#testing">test</A> program. </P>
<H2><A NAME="pointer"></A>5.11 Pointer arithmetic</H2>
<P CLASS="small"><A HREF="#err_hand">next</A> - <A HREF="#err_hand">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>Suppose you do something like </P>
<PRE>   int* y = new int[100];
   y += 200;          // y points to something outside the array
   // y is <B>never</B> accessed
</PRE>

<P>Then the standard says that the behaviour of the program is <I>undefined</I>
even if <TT>y</TT> is never accessed. (You are allowed to calculate a pointer
value one location beyond the end of the array). In practice, a program like
this does not cause any problems with any compiler I have come across and
no-one has reported any such problems to me. </P>
<P>However, this <I>error</I> is detected by Borland's <I>Code Guard</I>
bound's checker and this makes it very difficult to use this to use <I>Code
Guard</I> to detect other problems since the output is swamped by reports of
this <I>error</I>. </P>
<P>Now consider </P>
<PRE>   int* y = new int[100];
   y += 200;          // y points to something outside the array
   y -= 150;          // y points to something inside the array
   // y <B>is</B> accessed
</PRE>

<P>Again this is not strictly correct but does not seem to cause a problem. But
it is much more doubtful than the previous example. </P>
<P>I removed most instances of the second version of the problem from Newmat09.
Hopefully the remainder of these instances were removed from Newmat10. In addition, most instances of the first version of the
problem have also been fixed. </P>
<P>There is one exception. The interface to the <A HREF="#nric">Numerical
Recipes in C</A> does still contain the second version of the problem. This is
inevitable because of the way Numerical Recipes in C stores vectors and
matrices. If you are running the <A HREF="#testing">test program</A> with a
bounds checking program, edit <TT>tmt.h</TT> to disable the testing of the NRIC
interface. </P>
<P>The rule does does cause a problem for authors of matrix and
multidimensional array packages. If we want to run down a column of a matrix we
would like to do something like </P>
<PRE>   // set values of column 1
   Matrix A;
   ... set dimensions and put values in A
   Real* a = A.data();               // points to first element
   int nr = A.nrows();                // number of rows
   int nc = A.ncols();                // number of columns
   while (nr--)
   {
      *a = something to put in first element of row
      a += nc;                        // jump to next element of column
   }
</PRE>

<P>If the matrix has more than one column the last execution of <TT>a +=
nc;</TT> will run off the end of the space allocated to the matrix and we'll
get a bounds error report. </P>
<P>Instead we have to use a program like </P>
<PRE>   // set values of column 1
   Matrix A;
   ... set dimensions and put values in A
   Real* a = A.data();               // points to first element
   int nr = A.nrows();                // number of rows
   int nc = A.ncols();                // number of columns
   if (nr != 0)
   {
      for(;;)
      {
         *a = something to put in first element of row
         if (!(--nr)) break;
         a += nc;                     // jump to next element of column
      }
   }
</PRE>

<P>which is more complicated and consequently introduces more chance of error. 
</P>
<H2><A NAME="err_hand"></A>5.12 Error handling 
</H2>
<P CLASS="small"><A HREF="#sparse">next</A> - <A HREF="#sparse">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>The library now does have a moderately graceful exit from errors. One can
use either the simulated exceptions or the compiler supported exceptions. When
newmat08 was released (in 1995), compiler exception handling in the compilers I
had access to was unreliable. I recommended you used my simulated exceptions.
In 1997 compiler supported exceptions  seemed to work on a variety of
compilers - but not all compilers. This was still true in 2001. One compiler 
company was still having problems in 2003 (not sure about 2004). Try using the
compiler supported exceptions if you have a recent compiler, but if you are
getting strange crashes or errors try going back to my simulated exceptions.</P>
<P>The approach in the present library, attempting to simulate C++ exceptions,
is not completely satisfactory, but seems a good interim solution for those who
cannot use compiler supported exceptions. People who don't want exceptions in
any shape or form, can set the option to exit the program if an exception is
thrown. </P>
<P>The exception mechanism cannot clean-up objects explicitly created by new.
This must be explicitly carried out by the package writer or the package user.
I have not yet done this completely with the present package so occasionally a
little garbage may be left behind after an exception. I don't think this is a
big problem, but it is one that needs fixing. </P>
<H2><A NAME="sparse"></A>5.13 Sparse matrices</H2>
<P CLASS="small"><A HREF="#comp_mat">next</A> - <A HREF="#comp_mat">skip</A> -
<A HREF="#design">up</A> - <A HREF="#top">start</A></P>
<P>The library does not support sparse matrices. </P>
<P>For sparse matrices there is going to be some kind of structure vector. It
is going to have to be calculated for the results of expressions in much the
same way that types are calculated. In addition, a whole new set of row and
column operations would have to be written. </P>
<P>Sparse matrices are important for people solving large sets of differential
equations as well as being important for statistical and operational research
applications. </P>
<P>But there are packages being developed specifically for sparse matrices and
these might present the best approach, at least where sparse matrices are the
main interest. </P>
<H2><A NAME="comp_mat"></A>5.14 Complex matrices</H2>
<P CLASS="small"><a href="#function">next</a> - <a href="#function">skip</a> -
<A HREF="#design">up</A> - <A HREF="#top">start</A> </P>
<P>The package does not yet support matrices with complex elements. There are
at least two approaches to including these. One is to have matrices with
complex elements. </P>
<P>This probably means making new versions of the basic row and column
operations for Real*Complex, Complex*Complex, Complex*Real and similarly for
<TT>+</TT> and <TT>-</TT>. This would be OK, except that if I also want to do
this for sparse matrices, then when you put these together, the whole thing
will get out of hand. </P>
<P>The alternative is to represent a Complex matrix by a pair of Real matrices.
One probably needs another level of decoding expressions but I think it might
still be simpler than the first approach. But there is going to be a problem
with accessing elements and it does not seem possible to solve this in an
entirely satisfactory way. </P>
<P>Complex matrices are used extensively by electrical engineers and physicists
and really should be fully supported in a comprehensive package. </P>
<P>You can simulate most complex operations by representing <TT>Z = X + iY</TT>
by </P>
<PRE>    /  X   Y \
    \ -Y   X / 
</PRE>

<P>Most matrix operations will simulate the corresponding complex operation,
when applied to this matrix. But, of course, this matrix is essentially twice as big as you
would need with a genuine complex matrix library.</P>


<P>&nbsp;</P>


<h1><a name="function"></a>6. Function summary</h1>


<p class="small"><a href="#member_functions_1">next</a> - <a href="#changes">skip</a> -
<a href="#top">up</a> - <A HREF="#top">start</A> </p>


<p class="small"><a href="#member_functions_1">6.1 Member functions for matrices 
and matrix expressions</a><br>
<a href="#member_functions_2">6.2 Member functions for matrices</a><br>
<a href="#operators">6.3 Operators</a><br>
<a href="#global_functions_1">6.3 Global functions - newmat.h</a><br>
<a href="#global_functions_2">6.4 Global functions - newmatap.h</a><br>
<a href="#other_member">6.5 Other classes - member functions</a></p>


<p>This section lists member and global functions for matrices defined in <i>
newmat.h</i>. Where there are alternative names the lower-case non-capitalised 
versions are the preferred ones.</p>


<h2><a name="member_functions_1"></a>6.1 Member functions for matrices and matrix expressions</h2>
<p class="small"><a href="#member_functions_2">next</a> -
<a href="#member_functions_2">skip</a> - <a href="#function">up</a> - <A HREF="#top">start</A></p>
<p>Member functions for matrices and matrix expressions. These do not apply to 
<i>CroutMatrix</i> and <i>BandLUMatrix</i> 
unless explicitly noted.</p>


<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber1">
  <tr>
    <td width="33%"><b>Function group</b></td>
    <td width="33%"><b>function name</b></td>
    <td width="34%"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" rowspan="7"><a href="#unary">Unary operators</a><p>(see also
    <a href="#operators">operators</a>)</td>
    <td width="33%">.<b>t</b>()</td>
    <td width="34%">matrix transpose</td>
  </tr>
  <tr>
    <td width="33%">.<b>reverse</b>()<br>
    .Reverse()</td>
    <td width="34%">reverse order of elements (not band matrices)</td>
  </tr>
  <tr>
    <td width="33%">.<b>i</b>()</td>
    <td width="34%">invert matrix or solve (also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_rows</b>()</td>
    <td width="34%">sum elements in each row</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_columns</b>()</td>
    <td width="34%">sum  elements in each column</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_square_rows</b>()</td>
    <td width="34%">sum squares of elements in each row</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_square_columns</b>()</td>
    <td width="34%">sum squares of elements in each column</td>
  </tr>
  <tr>
    <td width="33%" rowspan="5"><a href="#ch_type">Change type</a></td>
    <td width="33%">.<b>as_row</b>()<br>
    .AsRow()</td>
    <td width="34%">interpret matrix body as a single row</td>
  </tr>
  <tr>
    <td width="33%">.<b>as_column</b>()<br>
    .AsColumn()</td>
    <td width="34%">interpret matrix body as a single column</td>
  </tr>
  <tr>
    <td width="33%">.<b>as_diagonal</b>()<br>
    .AsDiagonal()</td>
    <td width="34%">interpret matrix body as a diagonal matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>as_matrix</b>()<br>
    .AsMatrix()</td>
    <td width="34%">interpret matrix body as a rectangular matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>as_scalar</b>()<br>
    .AsScalar()</td>
    <td width="34%">convert 1x1 to Real</td>
  </tr>
  <tr>
    <td width="33%" rowspan="6"><a href="#submat">Submatrices</a></td>
    <td width="33%">.<b>submatrix</b>(int,int,int,int)<br>
    .SubMatrix(int,int,int,int)</td>
    <td width="34%">submatrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>sym_submatrix</b>(int,int)<br>
    .SymSubMatrix(int,int)</td>
    <td width="34%">submatrix with same row and column range</td>
  </tr>
  <tr>
    <td width="33%">.<b>row</b>(int)<br>
    .Row(int)</td>
    <td width="34%">select a row of a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>rows</b>(int,int)<br>
    .Rows(int,int)</td>
    <td width="34%">select a range of rows of a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>column</b>(int)<br>
    .Column(int,int)</td>
    <td width="34%">select a column of a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>columns</b>(int)<br>
    .Columns(int,int)</td>
    <td width="34%">select a range of columns of a matrix</td>
  </tr>
  <tr>
    <td width="33%" rowspan="12"><a href="#scalar2">Scalar functions - maxima &amp; 
    minima</a><p>(also global versions of maximum(), minimum(),&nbsp; 
    maximum_absolute_value(), mimimum_absolute_value())</td>
    <td width="33%">.<b>maximum_absolute_value</b>()<b><br>
    .</b>MaximumAbsoluteValue()</td>
    <td width="34%">maximum absolute value of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>maximum_absolute_value1</b>(int&amp;)<br>
    .MaximumAbsoluteValue1(int&amp;)</td>
    <td width="34%">maximum absolute value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>maximum_absolute_value2</b>(int&amp;,int&amp;)<br>
    .MaximumAbsoluteValue2(int&amp;,int&amp;)</td>
    <td width="34%">maximum absolute value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum_absolute_value</b>()<br>
    .MinimumAbsoluteValue()</td>
    <td width="34%">minimum absolute value of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum_absolute_value1</b>(int&amp;)<br>
    .MinimumAbsoluteValue1(int&amp;)</td>
    <td width="34%">minimum absolute value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum_absolute_value2</b>(int&amp;,int&amp;)<br>
    .MinimumAbsoluteValue2(int&amp;,int&amp;)</td>
    <td width="34%">minimum absolute value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>maximum</b>()<br>
    .Maximum()</td>
    <td width="34%">maximum value of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>maximum1</b>(int&amp;)<br>
    .Maximum1(int&amp;)</td>
    <td width="34%">maximum value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>maximum2</b>(int&amp;,int&amp;)<br>
    .Maximum2(int&amp;,int&amp;)</td>
    <td width="34%">maximum value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum</b>()<br>
    .Minimum()</td>
    <td width="34%">minimum value of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum1</b>(int&amp;)<br>
    .Minimum1(int&amp;)</td>
    <td width="34%">minimum value, return location</td>
  </tr>
  <tr>
    <td width="33%">.<b>minimum2</b>(int&amp;,int&amp;)<br>
    .Minimum2(int&amp;,int&amp;)</td>
    <td width="34%">minimum value, return location</td>
  </tr>
  <tr>
    <td width="33%" rowspan="9"><a href="#scalar3">Scalar functions - numerical</a><p>
    (also global versions of these functions)</td>
    <td width="33%">.<b>log_determinant</b>()<br>
    .LogDeterminant()</td>
    <td width="34%">natural logarithm of the determinant (also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%">.<b>determinant</b>()<br>
    .Determinant()</td>
    <td width="34%">determinant (also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_square</b>()<br>
    .SumSquare()</td>
    <td width="34%">sum of squares of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>norm_Frobenius</b>()<br>
    .norm_frobenius()<br>
    .NormFrobenius()</td>
    <td width="34%">square root of sum of squares of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum_absolute_value</b>()<br>
    .SumAbsoluteValue()</td>
    <td width="34%">sum of absolute values of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>sum</b>()<br>
    .Sum()</td>
    <td width="34%">sum of elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>trace</b>()<br>
    .Trace()</td>
    <td width="34%">trace of a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>norm1</b>()<br>
    .Norm1()</td>
    <td width="34%">maximum of sum of absolute values of elements of a column</td>
  </tr>
  <tr>
    <td width="33%">.<b>norm_infinity</b>()<br>
    .NormInfinity()</td>
    <td width="34%">maximum of sum of absolute values of elements of a row</td>
  </tr>
  <tr>
    <td width="33%"><a href="#scalar1">Scalar functions - size and shape</a></td>
    <td width="33%">.<b>bandwidth</b>()<br>
    .BandWidth()</td>
    <td width="34%">bandwidth of matrix (also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  </table>
<P>&nbsp;</P>
<h2><a name="member_functions_2"></a>6.2 Member functions for matrices</h2>
<p class="small"><a href="#operators">next</a> - <a href="#operators">skip</a> - <a href="#function">up</a> - <A HREF="#top">start</A></p>
<p>Member functions for matrices but not matrix expressions. These do not apply 
to CroutMatrix and BandLUMatrix unless explicitly noted.</p>
<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber2">
  <tr>
    <td width="33%"><b>Function group</b></td>
    <td width="33%"><b>function name</b></td>
    <td width="34%"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" rowspan="2"><a href="#elements">Element access</a><p>(see 
    also <a href="#operators">operators</a>)</td>
    <td width="33%">.<b>element</b>(int,int)</td>
    <td width="34%">access element - subscripts start at 0</td>
  </tr>
  <tr>
    <td width="33%">.<b>element</b>(int)</td>
    <td width="34%">access element - subscripts start at 0</td>
  </tr>
  <tr>
    <td width="33%" rowspan="2"><a href="#copy">Copying</a></td>
    <td width="33%">.<b>inject</b>(const GeneralMatrix&amp;)<br>
    .Inject(const GeneralMatrix&amp;)</td>
    <td width="34%">copy elements into a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>swap</b>(Matrix&amp;)</td>
    <td width="34%">swap bodies of two matrices of same type (also global 
	version, also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%" rowspan="9"><a href="#scalar1">Scalar functions - size and 
    shape</a><p>(also work with CroutMatrix and BandLUMatrix)</td>
    <td width="33%">.<b>type</b>()<br>
    .Type()</td>
    <td width="34%">type of a matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>nrows</b>()<br>
    .Nrows()</td>
    <td width="34%">number of rows</td>
  </tr>
  <tr>
    <td width="33%">.<b>ncols</b>()<br>
    .Ncols()</td>
    <td width="34%">number of columns</td>
  </tr>
  <tr>
    <td width="33%">.<b>size</b>()<br>
    .Storage()</td>
    <td width="34%">number of stored elements (including unused elements in band 
	matrices)</td>
  </tr>
  <tr>
    <td width="33%">.<b>size2</b>()</td>
    <td width="34%">size of second array (BandLUMatrix only)</td>
  </tr>
  <tr>
    <td width="33%">.<b>data</b>()<br>
    .Store()</td>
    <td width="34%">pointer to stored elements</td>
  </tr>
  <tr>
    <td width="33%">.<b>const_data</b>()</td>
    <td width="34%">constant pointer to stored elements</td>
  </tr>
  <tr>
    <td width="33%"><b>.const_data2</b>()</td>
    <td width="34%">constant pointer to second array (BandLUMatrix only)</td>
  </tr>
  <tr>
    <td width="33%"><b>.const_data_indx</b>()</td>
    <td width="34%">constant pointer to row swap array (CroutMatrix and BandLUMatrix 
	only)</td>
  </tr>
  <tr>
    <td width="33%" rowspan="3"><a href="#scalar3">Scalar functions - numerical</a></td>
    <td width="33%">.<b>is_zero</b>()<br>
    .IsZero()</td>
    <td width="34%">test all elements are exactly zero (also global version)</td>
  </tr>
  <tr>
    <td width="33%">.<b>is_singular</b>()<br>
    .IsSingular()</td>
    <td width="34%">test for exact singularity (CroutMatrix and BandLUMatrix 
	only)</td>
  </tr>
  <tr>
    <td width="33%"><b>.even_exchanges</b>()</td>
    <td width="34%">true if there have been an even number of row exchanges (CroutMatrix and BandLUMatrix 
	only)</td>
  </tr>
  <tr>
    <td width="33%" rowspan="4"><a href="#memory">Memory management</a></td>
    <td width="33%">.<b>release</b>()<br>
    .Release()</td>
    <td width="34%">release memory after next operation</td>
  </tr>
  <tr>
    <td width="33%">.<b>release</b>(int)<br>
    .Release(int)</td>
    <td width="34%">release memory after specified number of operations</td>
  </tr>
  <tr>
    <td width="33%">.<b>release_and_delete</b>()<br>
    .ReleaseAndDelete()</td>
    <td width="34%">delete after next operation</td>
  </tr>
  <tr>
    <td width="33%">.<b>for_return</b>()<br>
    .ForReturn()</td>
    <td width="34%"><i>place in an envelope</i> for efficient return from a function</td>
  </tr>
  <tr>
    <td width="33%" rowspan="7"><a href="#dimen">Change dimensions</a></td>
    <td width="33%">.<b>resize</b>(int)<br>
    .ReSize(int)</td>
    <td width="34%">change the dimensions (vectors and square matrices)</td>
  </tr>
  <tr>
    <td width="33%">.<b>resize</b>(int,int)<br>
    .ReSize(int,int)</td>
    <td width="34%">change the dimensions (non-square matrices, triangular band 
    matrices and symmetric band matrices)</td>
  </tr>
  <tr>
    <td width="33%">.<b>resize</b>(int,int,int)<br>
    .ReSize(int,int,int)</td>
    <td width="34%">change the dimensions (band matrices)</td>
  </tr>
  <tr>
    <td width="33%">.<b>resize</b>(const GeneralMatrix&amp;)<br>
    .ReSize(const GeneralMatrix&amp;)</td>
    <td width="34%">change dimensions to match those of another matrix</td>
  </tr>
  <tr>
    <td width="33%">.<b>cleanup</b>()<br>
    .CleanUp()</td>
    <td width="34%">resize to 0x0</td>
  </tr>
  <tr>
    <td width="33%"><b>.resize_keep</b>(int)</td>
    <td width="34%">change the dimensions, keep values (vectors and square matrices, 
    not band)</td>
  </tr>
  <tr>
    <td width="33%">.<b>resize_keep</b>(int,int)</td>
    <td width="34%">change the dimensions , keep values (non-square matrices)</td>
  </tr>
</table>
<P>&nbsp;</P>
<h2><a name="operators"></a>6.3 Operators</h2>
<p class="small"><a href="#global_functions_1">next</a> -
<a href="#global_functions_1">skip</a> - <a href="#function">up</a> - <A HREF="#top">start</A></p>
<P>Operators for matrices and matrix expressions. These do not apply to CroutMatrix and BandLUMatrix 
unless explicitly noted.</P>


<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber1">
  <tr>
    <td width="33%"><b>Function group</b></td>
    <td width="33%"><b>function name</b></td>
    <td width="34%"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" rowspan="3"><a href="#elements">Element access</a><p>
    (matrices only, not functions of a matrix)</td>
    <td width="33%"><tt>()</tt></td>
    <td width="34%">access element - subscripts start at 1</td>
  </tr>
  <tr>
    <td width="33%"><tt>()</tt></td>
    <td width="34%">access element - subscripts start at 1</td>
  </tr>
  <tr>
    <td width="33%"><tt>[]</tt></td>
    <td width="34%">access element C style - subscripts start at zero; if <i>
    SETUP_C_SUBSCRIPTS</i> is defined.</td>
  </tr>
  <tr>
    <td width="33%"><a href="#unary">Unary operators</a></td>
    <td width="33%"><tt>-</tt></td>
    <td width="34%">change sign of elements</td>
  </tr>
  <tr>
    <td width="33%" rowspan="7"><a href="#binary">Binary operators</a></td>
    <td width="33%"><tt>+</tt>, <tt>+=</tt></td>
    <td width="34%">add matrices</td>
  </tr>
  <tr>
    <td width="33%"><tt>-</tt>, <tt>-=</tt></td>
    <td width="34%">subtract matrices</td>
  </tr>
  <tr>
    <td width="33%"><tt>*</tt>, <tt>*=</tt></td>
    <td width="34%">matrix multiplication</td>
  </tr>
  <tr>
    <td width="33%"><tt>|</tt>, <tt>|=</tt></td>
    <td width="34%">horizontal concatenation</td>
  </tr>
  <tr>
    <td width="33%"><tt>&amp;</tt>, <tt>&amp;=</tt></td>
    <td width="34%">vertical stacking</td>
  </tr>
  <tr>
    <td width="33%"><tt>==</tt></td>
    <td width="34%">test for <i>exact equality </i>(also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%"><tt>!=</tt></td>
    <td width="34%">test for inequality (i.e. not <i>exact equality</i>, also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%" rowspan="4"><a href="#matscal">Matrix and scalar</a></td>
    <td width="33%"><tt>+</tt>, <tt>+=</tt></td>
    <td width="34%">add <i>Real</i> to matrix</td>
  </tr>
  <tr>
    <td width="33%"><tt>-</tt>, <tt>-=</tt></td>
    <td width="34%">subtract <i>Real</i> from matrix; subtract matrix from <i>Real</i></td>
  </tr>
  <tr>
    <td width="33%"><tt>*</tt>, <tt>*=</tt></td>
    <td width="34%">multiply matrix by <i>Real</i></td>
  </tr>
  <tr>
    <td width="33%"><tt>/</tt>, <tt>/=</tt></td>
    <td width="34%">divide matrix by <i>Real</i></td>
  </tr>
  <tr>
    <td width="33%" rowspan="3"><a href="#copy">Copying</a></td>
    <td width="33%"><tt>=</tt></td>
    <td width="34%">copy matrix (error if there is loss of data, <i>&nbsp;</i>also works with CroutMatrix and BandLUMatrix)</td>
  </tr>
  <tr>
    <td width="33%"><tt>=</tt></td>
    <td width="34%">copy <i>Real</i> 
    to all elements</td>
  </tr>
  <tr>
    <td width="33%"><tt>&lt;&lt;</tt></td>
    <td width="34%">copy matrix (no error if there is loss of data)</td>
  </tr>
  <tr>
    <td width="33%"><a href="#entering">Enter values</a></td>
    <td width="33%"><tt>&lt;&lt;</tt></td>
    <td width="34%">enter list of values into matrix</td>
  </tr>
  <tr>
    <td width="33%"><a href="#output">Output</a><p>(header in <i>newmatio.h</i>)</td>
    <td width="33%"><tt>&lt;&lt;</tt></td>
    <td width="34%">print matrix to file</td>
  </tr>
  </table>
<P>&nbsp;</P>

<h2><a name="global_functions_1"></a>6.4 Global functions - newmat.h</h2>
<p class="small"><a href="#global_functions_2">next</a> -
<a href="#global_functions_2">skip</a> -
<a href="#function">up</a> - <A HREF="#top">start</A></p>
<p>Operators for matrices and matrix expressions. These do not apply to CroutMatrix and BandLUMatrix.</p>
<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber3">
  <tr>
    <td width="33%"><b>Function group</b></td>
    <td width="33%"><b>function name</b></td>
    <td width="34%"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" rowspan="5"><a href="#binary">Binary operators</a><p>(see 
    also <a href="#operators">operators</a>)</td>
    <td width="33%"><b>SP</b>(const BaseMatrix&amp;, const BaseMatrix&amp;)</td>
    <td width="34%">element-wise product of two matrices</td>
  </tr>
  <tr>
    <td width="33%"><b>KP</b>(const BaseMatrix&amp;, const BaseMatrix&amp;)</td>
    <td width="34%">Kronecker product of two matrices</td>
  </tr>
  <tr>
    <td width="33%"><b>crossproduct</b>(const Matrix&amp;, const Matrix&amp;)<br>
    CrossProduct(const Matrix&amp;, const Matrix&amp;)</td>
    <td width="34%">cross product of two 3x1 or 1x3 matrices or vectors.</td>
  </tr>
  <tr>
    <td width="33%"><b>crossproduct_rows</b>(const Matrix&amp;, const Matrix&amp;)<br>
    CrossProductRows(const Matrix&amp;, const Matrix&amp;)</td>
    <td width="34%">row-wise cross product</td>
  </tr>
  <tr>
    <td width="33%"><b>crossproduct_columns</b>(const Matrix&amp;, const Matrix&amp;)<br>
    CrossProductColumns(const Matrix&amp;, const Matrix&amp;)</td>
    <td width="34%">column-wise cross product</td>
  </tr>
  <tr>
    <td width="33%"><a href="#scalar3">Scalar functions - numerical</a></td>
    <td width="33%"><b>dotproduct</b>(const Matrix&amp;, const Matrix&amp;)<b><br>
    </b>DotProduct(const Matrix&amp;, const Matrix&amp;)</td>
    <td width="34%">dot product of two vectors</td>
  </tr>
  </table>

<P>&nbsp;</P>


<h2><a name="global_functions_2"></a>6.5 Global functions - newmatap.h</h2>
<p class="small"><a href="#other_member">next</a> - <a href="#other_member">skip</a> -
<a href="#function">up</a> - <A HREF="#top">start</A></p>

<P>Advanced operators for matrices and matrix expressions. These do not apply to CroutMatrix and BandLUMatrix.</P>


<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber1">
  <tr>
    <td width="33%" height="19"><b>Function group</b></td>
    <td width="33%" height="19"><b>function name</b></td>
    <td width="34%" height="19"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" rowspan="11" height="195"><a href="#qr">QR transform</a></td>
    <td width="33%" height="38"><b>QRZT</b>(Matrix&amp;, LowerTriangularMatrix&amp;)</td>
    <td width="34%" height="38">transposed version of QRZ transform</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>QRZT</b>(const Matrix&amp;, Matrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">transposed version of QRZ transform - solve part</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>QRZT</b>(Matrix&amp;, Matrix&amp;, LowerTriangularMatrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="19">both of previous two lines</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>QRZ</b>(Matrix&amp;, UpperTriangularMatrix&amp;)</td>
    <td width="34%" height="19">QRZ transform</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>QRZ</b>(const Matrix&amp;, Matrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="19">QRZ transform - solve part</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>QRZ</b>(Matrix&amp;, Matrix&amp;, UpperTriangularMatrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">both of previous two lines</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>updateQRZT</b>(Matrix&amp;, LowerTriangularMatrix&amp;)<br>
    UpdateQRZT(Matrix&amp;, LowerTriangularMatrix&amp;)</td>
    <td width="34%" height="38">add extra columns to transposed version of QRZ 
    transform</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>updateQRZ</b>(Matrix&amp;, UpperTriangularMatrix&amp;)<br>
    UpdateQRZ(Matrix&amp;, UpperTriangularMatrix&amp;)</td>
    <td width="34%" height="38">add extra rows to QRZ transform</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>updateQRZ</b>(const Matrix&amp;, Matrix&amp;, Matrix&amp;)<br>
	UpdateQRZ(const Matrix&amp;, Matrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">combine the results of the solve parts of QRZ 
	transforms on two blocks of data</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>updateQRZ</b>(UpperTriangularMatrix&amp;, UpperTriangularMatrix&amp;)<br>
	UpdateQRZ(UpperTriangularMatrix&amp;, UpperTriangularMatrix&amp;)</td>
    <td width="34%" height="38">combine the results of QRZ transforms on two 
	blocks of data</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>updateQRZ</b>(const UpperTriangularMatrix&amp;, Matrix&amp;, Matrix&amp;)<br>
	UpdateQRZ(const UpperTriangularMatrix&amp;, Matrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">combine the results of the solve parts of QRZ 
	transforms on two blocks of data</td>
  </tr>
  <tr>
    <td width="33%" height="38"><a href="#extend">Extend orthonormal set</a></td>
    <td width="33%" height="38"><b>extend_orthonormal</b>(Matrix&amp;, int)</td>
    <td width="34%" height="38">extend a set of orthonormal columns</td>
  </tr>
  <tr>
    <td width="33%" rowspan="2" height="77"><a href="#cholesky">Cholesky 
    decomposition</a></td>
    <td width="33%" height="38"><b>Cholesky</b>(const SymmetricMatrix&amp;)</td>
    <td width="34%" height="38">Cholesky decomposition of symmetric matrix</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>Cholesky</b>(const SymmetricBandMatrix&amp;)</td>
    <td width="34%" height="38">Cholesky decomposition of symmetric band matrix</td>
  </tr>
  <tr>
    <td width="33%" rowspan="4" height="155"><a href="#upd_chol">Update Cholesky 
    decomposition</a></td>
    <td width="33%" height="38"><b>update_Cholesky </b>(UpperTriangularMatrix&amp;, 
    RowVector)<br>
    UpdateCholesky<b> </b>(UpperTriangularMatrix&amp;, 
    RowVector)</td>
    <td width="34%" height="38">add extra row to Cholesky/QR decomposition</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>downdate_Cholesky </b>(UpperTriangularMatrix&amp;, 
    RowVector)<br>
    DowndateCholesky<b> </b>(UpperTriangularMatrix&amp;, 
    RowVector)</td>
    <td width="34%" height="38">remove row from Cholesky/QR decomposition</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>right_circular_update_Cholesky </b>(UpperTriangularMatrix&amp;, 
    int, int)<br>
    RightCircularUpdateCholesky<b> </b>(UpperTriangularMatrix&amp;, 
    int, int)</td>
    <td width="34%" height="38">rearrange columns in Cholesky/QR decomposition</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>left_circular_update_Cholesky </b>(UpperTriangularMatrix&amp;, 
    int, int)<br>
    LeftCircularUpdateCholesky<b> </b>(UpperTriangularMatrix&amp;, 
    int, int)</td>
    <td width="34%" height="38">rearrange columns in Cholesky/QR decomposition</td>
  </tr>
  <tr>
    <td width="33%" rowspan="3" height="99"><a href="#svd">Singular value 
    decomposition</a></td>
    <td width="33%" height="38"><b>SVD</b>(const Matrix&amp;, DiagonalMatrix&amp;, 
    Matrix&amp;, Matrix&amp;, bool, bool)</td>
    <td width="34%" height="38">singular value decomposition - get U and V</td>
  </tr>
  <tr>
    <td width="33%" height="21"><b>SVD</b>(const Matrix&amp;, DiagonalMatrix&amp;)</td>
    <td width="34%" height="21">SVD decomposition - get just singular values</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>SVD</b>(const Matrix&amp; A, DiagonalMatrix&amp; D, 
    Matrix&amp;, bool)</td>
    <td width="34%" height="38">SVD decomposition - get U</td>
  </tr>
  <tr>
    <td width="33%" rowspan="7" height="234"><a href="#evalues">Eigenvalue 
    decomposition of a symmetric matrix</a></td>
    <td width="33%" height="19"><b>Jacobi</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;)</td>
    <td width="34%" height="19">Jacobi eigenvalue decomposition - get only eigenvalues</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>Jacobi</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, SymmetricMatrix&amp;)</td>
    <td width="34%" height="38">Jacobi eigenvalue decomposition - get only 
    eigenvalues</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>Jacobi</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">Jacobi eigenvalue decomposition - get 
    eigenvalues and eigenvectors</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>Jacobi</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, SymmetricMatrix&amp;, Matrix&amp;, bool)</td>
    <td width="34%" height="38">Jacobi eigenvalue decomposition - get 
    eigenvalues and eigenvectors</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>eigenvalues</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;)<br>
    EigenValues(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;)</td>
    <td width="34%" height="19">Householder eigenvalue decomposition - get only 
    eigenvalues </td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>eigenvalues</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, SymmetricMatrix&amp;)<br>
    EigenValues(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, SymmetricMatrix&amp;)</td>
    <td width="34%" height="38">Householder eigenvalue decomposition with back 
    transform - get only eigenvalues</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>eigenvalues</b>(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, Matrix&amp;)<br>
    EigenValues(const SymmetricMatrix&amp;, 
    DiagonalMatrix&amp;, Matrix&amp;)</td>
    <td width="34%" height="38">Householder eigenvalue decomposition - get 
    eigenvalues and eigenvectors</td>
  </tr>
  <tr>
    <td width="33%" rowspan="2" height="39"><a href="#sorting">Sorting</a></td>
    <td width="33%" height="19"><b>sort_ascending</b>(GeneralMatrix&amp;)<br>
    SortAscending(GeneralMatrix&amp;)</td>
    <td width="34%" height="19">ascending sort</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>sort_descending</b>(GeneralMatrix&amp;)<br>
    SortDescending(GeneralMatrix&amp;)</td>
    <td width="34%" height="19">descending sort</td>
  </tr>
  <tr>
    <td width="33%" rowspan="6" height="233"><a href="#fft">Fast Fourier 
    transform</a></td>
    <td width="33%" height="38"><b>FFT</b>(const ColumnVector&amp;, const 
    ColumnVector&amp;, ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="38">fast Fourier transform</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>FFTI</b>(const ColumnVector&amp;, const 
    ColumnVector&amp;, ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="38">fast Fourier transform - inverse</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>RealFFT</b>(const ColumnVector&amp;, ColumnVector&amp;, 
    ColumnVector&amp;)</td>
    <td width="34%" height="38">FFT of real vector</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>RealFFTI</b>(const ColumnVector&amp;, const 
    ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="38">FFT of real vector - inverse</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>FFT2</b>(const Matrix&amp; U, const Matrix&amp; V, 
    Matrix&amp; X, Matrix&amp; Y)</td>
    <td width="34%" height="38">two dimensional FFT</td>
  </tr>
  <tr>
    <td width="33%" height="38"><b>FFT2I</b>(const Matrix&amp; U, const Matrix&amp; V, 
    Matrix&amp; X, Matrix&amp; Y)</td>
    <td width="34%" height="38">two dimensional FFT - inverse</td>
  </tr>
  <tr>
    <td width="33%" rowspan="8" height="80"><a href="#trigtran">Fast 
    trigonometric transform</a></td>
    <td width="33%" height="19"><b>DCT_II</b>(const ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="19">type II discrete cosine transform</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DCT_II_inverse</b>(const ColumnVector&amp;, 
    ColumnVector&amp;)</td>
    <td width="34%" height="19">type II discrete cosine transform - inverse</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DST_II</b>(const ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="19">type II discrete sine transform</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DST_II_inverse</b>(const ColumnVector&amp;, 
    ColumnVector&amp;)</td>
    <td width="34%" height="19">type II discrete sine transform - inverse</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DCT</b>(const ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="19">discrete cosine transform</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DCT_inverse</b>(const ColumnVector&amp;, 
    ColumnVector&amp;)</td>
    <td width="34%" height="19">discrete cosine transform - inverse</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DST</b>(const ColumnVector&amp;, ColumnVector&amp;)</td>
    <td width="34%" height="19">discrete sine transform</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>DST_inverse</b>(const ColumnVector&amp;, 
    ColumnVector&amp;)</td>
    <td width="34%" height="19">discrete sine transform - inverse</td>
  </tr>
  <tr>
    <td width="33%" rowspan="5" height="79">
	<a href="#misc_fn">Helmert transform</a></td>
    <td width="33%" height="19"><b>Helmert</b>(int, bool=false)</td>
    <td width="34%" height="19">return Helmert transform matrix</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>Helmert</b>(const ColumnVector&amp;, bool=false)<br>
	<b>Helmert</b>(const Matrix&amp;, bool=false)</td>
    <td width="34%" height="19">multiply by Helmert transform matrix</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>Helmert</b>(int, int, bool=false)</td>
    <td width="34%" height="19">return column of Helmert transform matrix</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>Helmert_transpose</b>(const ColumnVector&amp;, bool=false)<br>
	<b>Helmert_transpose</b>(const Matrix&amp;, bool=false)</td>
    <td width="34%" height="19">multiply by transpose of Helmert transform 
	matrix</td>
  </tr>
  <tr>
    <td width="33%" height="19"><b>Helmert_transpose</b>(const ColumnVector&amp;, int, bool=false)</td>
    <td width="34%" height="19">multiply by transpose of Helmert transform 
	matrix, return one element of result</td>
  </tr>
  </table>

<P>&nbsp;</P>
<h2><a name="other_member"></a>6.6 Other classes - member functions</h2>
<p class="small"><a href="#changes">next</a> - <a href="#changes">skip</a> -
<a href="#function">up</a> - <A HREF="#top">start</A></p>

<P>&nbsp;</P>


<table border="1" cellpadding="0" cellspacing="1" style="border-collapse: collapse" bordercolor="#111111" width="100%" id="AutoNumber1">
  <tr>
    <td width="33%" height="19"><b>Class</b></td>
    <td width="33%" height="19"><b>function name</b></td>
    <td width="34%" height="19"><b>description</b></td>
  </tr>
  <tr>
    <td width="33%" height="95" rowspan="5"><a href="#scalar3">LogAndSign</a></td>
    <td width="33%" height="19">.<b>pow_eq</b>(int)<br>
    .PowEq(int)</td>
    <td width="34%" height="19">raise to a power</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>change_sign</b>()<br>
    .ChangeSign()</td>
    <td width="34%" height="19">change sign</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>log_value</b>()<br>
    .LogValue()</td>
    <td width="34%" height="19">return the natural logarithm of the value</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>sign</b>()<br>
    .Sign()</td>
    <td width="34%" height="19">return the sign</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>value</b>()<br>
    .Value()</td>
    <td width="34%" height="19">return the value (no log transform)</td>
  </tr>
  <tr>
    <td width="33%" height="76" rowspan="4"><a href="#scalar1">MatrixType</a></td>
    <td width="33%" height="19">.<b>value</b>()<br>
    .Value()</td>
    <td width="34%" height="19">return the value (as character string)</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>is_diagonal</b>()</td>
    <td width="34%" height="19">has diagonal attribute</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>is_symmetric</b>() </td>
    <td width="34%" height="19">has symmetric attribute</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>is_band</b>()</td>
    <td width="34%" height="19">has band attribute</td>
  </tr>
  <tr>
    <td width="33%" height="38" rowspan="2"><a href="#scalar1">MatrixBandWidth</a></td>
    <td width="33%" height="19">.<b>upper</b>()<br>
    .Upper()</td>
    <td width="34%" height="19">return upper band width</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>lower</b>()<br>
    .Lower()</td>
    <td width="34%" height="19">return lower bandwidth</td>
  </tr>
  <tr>
    <td width="33%" height="114" rowspan="6"><a href="#SimpleIntArray">SimpleIntArray</a></td>
    <td width="33%" height="19">.<b>size</b>()<br>
    .Size()</td>
    <td width="34%" height="19">return size of array</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>data</b>()<br>
    .Data()</td>
    <td width="34%" height="19">return a pointer to the data</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>const_data</b>()</td>
    <td width="34%" height="19">return a constant pointer to the data</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>resize</b>()<br>
    .Resize()</td>
    <td width="34%" height="19">change the size of an array</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>resize_keep</b>()<br>
    .resize(true)<br>
    .Resize(true)</td>
    <td width="34%" height="19">change the size of an array - keep the values</td>
  </tr>
  <tr>
    <td width="33%" height="19">.<b>cleanup</b>()<br>
    .CleanUp()</td>
    <td width="34%" height="19">resize to zero</td>
  </tr>
</table>


<P>&nbsp;</P>


<h1><A NAME="changes"></A>7. Change history</h1>
<P CLASS="small"><a href="#problem">next</a> - <a href="#problem">skip</a> -
<a href="#top">up</a> - <a href="#top">start</a></P>
<P><b>Newmat11 - November, 2008:</b></P>
<P CLASS="small">Remove work-arounds for older compilers, Borland Builder 6 and 
Open Watcom compatibility, SquareMatrix, 
load from array of ints, crossproducts, Cholesky and QRZ update functions, swap 
functions, FFT2, access to arrays in traditional C functions, SimpleIntArray 
class, compatibility with Numerical Recipes in C++, sum_rows(), sum_columns(), sum_squares_rows() and 
sum_squares_columns() functions, extend_orthogonal function, resize_keep 
function, speed-ups and 
bug-fixes,  change to lower case for functions, can copy CroutMatrix, 
BandLUMatrix, Helmert transform, compatibility fix for G++ 4.1, start inserting 
comments for Doxygen, SP_eq, scientific format for output, fix for 64 bits.</P>
<P><b>Newmat10A - October, 2002, Newmat10B - January 2005:</b></P>
<P CLASS="small">Fix error in Kronecker product; fixes for Intel and GCC3 
compilers.</P>
<P><B>Newmat10 - January, 2002:</B> </P>
<P CLASS="small">Improve compatibility with GCC, fix errors in FFT and
GenericMatrix, update simulated exceptions, maxima, minima, determinant, dot
product and Frobenius norm functions, update make files for CC and GCC, faster
FFT, <TT>A.ReSize(B)</TT>, fix pointer arithmetic, <TT>&lt;&lt;</TT> for loading 
data into rows, IdentityMatrix, Kronecker product, sort singular values.</P>
<P><B>Newmat09 - September, 1997:</B> </P>
<P CLASS="small">Operator <TT>==</TT>, <TT>!=</TT>, <TT>+=</TT>, <TT>-=</TT>,
<TT>*=</TT>, <TT>/=</TT>, <TT>|=</TT>, <TT>&amp;=</TT>. Follow new rules for
<I>for (int i; ... )</I> construct. Change Boolean, TRUE, FALSE to bool, true,
false. Change ReDimension to ReSize. SubMatrix allows zero rows and columns.
Scalar <TT>+</TT>, <TT>-</TT> or <TT>*</TT> matrix is OK. Simplify simulated
exceptions. Fix non-linear programs for AT&amp;T compilers. Dummy inequality
operators. Improve internal row/column operations. Improve matrix LU
decomposition. Improve sort. Reverse function. IsSingular function. Fast trig
transforms. Namespace definitions. </P>
<P><B>Newmat08A - July, 1995:</B> </P>
<P CLASS="small">Fix error in SVD. </P>
<P><B>Newmat08 - January, 1995:</B> </P>
<P CLASS="small">Corrections to improve compatibility with Gnu, Watcom.
Concatenation of matrices. Elementwise products. Option to use compilers
supporting exceptions. Correction to exception module to allow global
declarations of matrices. Fix problem with inverse of symmetric matrices. Fix
divide-by-zero problem in SVD. Include new QR routines. Sum function.
Non-linear optimisation. GenericMatrices. </P>
<P><B>Newmat07 - January, 1993</B> </P>
<P CLASS="small">Minor corrections to improve compatibility with Zortech,
Microsoft and Gnu. Correction to exception module. Additional FFT functions.
Some minor increases in efficiency. Submatrices can now be used on RHS of =.
Option for allowing C type subscripts. Method for loading short lists of
numbers. </P>
<P><B>Newmat06 - December 1992:</B> </P>
<P CLASS="small">Added band matrices; 'real' changed to 'Real' (to avoid
potential conflict in complex class); Inject doesn't check for no loss of
information; fixes for AT&amp;T C++ version 3.0; real(A) becomes A.AsScalar();
CopyToMatrix becomes AsMatrix, etc; .c() is no longer required (to be deleted
in next version); option for version 2.1 or later. Suffix for include files
changed to .h; BOOL changed to Boolean (BOOL doesn't work in g++ v 2.0);
modifications to allow for compilers that destroy temporaries very quickly;
(Gnu users - see the section of compilers). Added CleanUp,
LinearEquationSolver, primitive version of exceptions. </P>
<P><B>Newmat05 - June 1992:</B> </P>
<P CLASS="small">For private release only </P>
<P><B>Newmat04 - December 1991:</B> </P>
<P CLASS="small">Fix problem with G++1.40, some extra documentation  
</P>
<P><B>Newmat03 - November 1991:</B> </P>
<P CLASS="small">Col and Cols become Column and Columns. Added Sort, SVD,
Jacobi, Eigenvalues, FFT, real conversion of 1x1 matrix, <I>Numerical Recipes
in C</I> interface, output operations, various scalar functions. Improved
return from functions. Reorganised setting options in &quot;include.hxx&quot;.
</P>
<P><B>Newmat02 - July 1991:</B> </P>
<P CLASS="small">Version with matrix row/column operations and numerous
additional functions. </P>
<P><B>Matrix - October 1990:</B> </P>
<P CLASS="small">Early version of package. </P>

<P CLASS="small">&nbsp;</P>
<H2><A NAME="problem"></A>8. Problem report form
</H2>

<P CLASS="small"><a href="#top">next</a> - <a href="#top">skip</a> - <a href="#top">up</a> -
<a href="#top">start</a></P>

<P>Copy and paste this to your editor; fill it out and email to robert <b>at</b> statsresearch.co.nz </P>
<P>But first look in my web page <a href="http://www.robertnz.net">
http://www.robertnz.net</a> to see if the bug has
already been reported. </P>
<PRE> Version: ............... newmat11 (20 November 2008)
 Your email address: ....
 Today's date: ..........
 Your machine: ..........
 Operating system: ......
 Compiler &amp; version: ....
 Compiler options
   (eg GUI or console)...
 Describe the problem - attach examples if possible:



</PRE>

<PRE>-----------------------------------------------------------
</PRE>

<P CLASS="small">&nbsp;</P>

<P CLASS="small">&nbsp;</P>
<P CLASS="small"><A HREF="#top">return to top</A><BR>
<A HREF="ol_doc.htm">return to online documentation page</A></P>
<P>&nbsp;</P>
<P>&nbsp;</P>
<P>&nbsp;</P>
<P>&nbsp;</P>
<P>&nbsp;</P>
<P>&nbsp;</P>
<P>&nbsp;</P>
</BODY>
</HTML>